On the theory of noncovalent binding

On Computing Equilibrium Binding Constants for Protein-Protein Association in Membranes

Protein association in lipid membranes is fundamental to membrane protein function and of great biomedical relevance. All-atom and coarse-grained models have been extensively used to understand the protein-protein interactions in the membrane and to compute equilibrium association constants. However, slow translational and rotational diffusion of protein in membrane presents challenges to the effective sampling of conformations defining the ensembles of free and bound states contributing to the association equilibrium and the free energy of dimerization. We revisit the homodimerization equilibrium of the TM region of glycophorin A. Conformational sampling is performed using umbrella sampling along previously proposed one-dimensional collective variables and compared with sampling over a two-dimensional collective variable space using the MARTINI v2.2 force field. We demonstrate that the one-dimensional collective variables suffer from restricted sampling of the native homodimer conformations leading to a biased free energy landscape. Conversely, simulations along the two-dimensional collective variable effectively characterize the thermodynamically relevant native and non-native interactions contributing to the association equilibrium. These results demonstrate the challenges associated with accurately characterizing binding equilibria when multiple poses contribute to the bound state ensemble.

Free Energy-Based Computational Methods for the Study of Protein-Peptide Binding Equilibria

This chapter discusses the theory and application of physics-based free energy methods to estimate protein-peptide binding free energies. It presents a statistical mechanics formulation of molecular binding, which is then specialized in three methodologies: (1) alchemical absolute binding free energy estimation with implicit solvation, (2) alchemical relative binding free energy estimation with explicit solvation, and (3) potential of mean force binding free energy estimation. Case studies of protein-peptide binding application taken from the recent literature are discussed for each method.

Building intuition for binding free energy calculations: Bound state definition, restraints, and symmetry

The theory behind computation of absolute binding free energies using explicit-solvent molecular simulations is well-established, yet somewhat complex, with counter-intuitive aspects. This leads to frequent frustration, common misconceptions, and sometimes erroneous numerical treatment. To improve this, we present the main practically relevant segments of the theory with constant reference to physical intuition. We pinpoint the role of the implicit or explicit definition of the bound state (or the binding site) to make a robust link between an experimental measurement and a computational result. We clarify the role of symmetry and discuss cases where symmetry number corrections have been misinterpreted. In particular, we argue that symmetry corrections as classically presented are a source of confusion and could be advantageously replaced by restraint free energy contributions. We establish that contrary to a common intuition, partial or missing sampling of some modes of symmetric bound states does not affect the calculated decoupling free energies. Finally, we review these questions and pitfalls in the context of a few common practical situations: binding to a symmetric receptor (equivalent binding sites), binding of a symmetric ligand (equivalent poses), and formation of a symmetric complex, in the case of homodimerization.

Noncovalent Association and Partitioning of Some Perfume Components at Infinite Dilution with Myelin Basic Protein Pseudophase in Normal Saline

Myelin basic protein (MBP), one of the major protein constituents of the myelin sheath, possesses unique ligand-binding features. We present a novel equilibrium headspace gas chromatographic technique to examine the thermodynamics of noncovalent interactions between common perfume components: Lilial, Hedione, Hexylcinnamic aldehyde, and Versalide with MBP monomers and its hexameric MBP-pseudophase. A general theoretical model is used to calculate the critical aggregation concentration (cac) of MBP, perfume component binding constants with monomeric MBP, K₁₁, and MBP as pseudophase, K_n1, and free energies for perfume component binding with monomeric MBP, ΔG_b,11, and MBP as pseudophase, ΔG_b,n1. In addition, the pseudophase-water partition coefficients, K_x, the free energies of transfer of perfume from bulk water to the MBP-pseudophase, ΔG_t, and the intra-aggregate activity coefficients, γ_m^∞, at infinite dilution were also determined. The cac value measured by the method of fractional distribution is a unique and precise approach in understanding the aggregation phenomenon. Within the experimental error, the 1:1 binding free energies did not differ by more than 1 kJ/mol among the perfume components but favored the MBP pseudophase binding by 6 kJ/mol. Therefore, that protein aggregation can enhance the binding of small molecules is probably a general conclusion. While the magnitudes of K₁₁, K_n1, ΔG_b, K_x, and ΔG_t show weak trends, the γ_m^∞ values show a strong and distinct trend in interaction, spanning 4 orders of magnitude among the perfume components.

Methodological uncertainties in drug-receptor binding free energy predictions based on classical molecular dynamics

Computational approaches are becoming an essential tool in modern drug design and discovery, with fast compound triaging using a combination of machine learning and docking techniques followed by molecular dynamics binding free energies assessment using alchemical techniques. The traditional MD-based alchemical free energy perturbation (FEP) method faces severe sampling issues that may limits its reliability in automated workflows. Here we review the major sources of uncertainty in FEP protocols for drug discovery, showing how the sampling problem can be effectively tackled by switching to nonequilibrium alchemical techniques.

Alchemical free energy calculations via metadynamics: Application to the theophylline-RNA aptamer complex

We propose a computational workflow for robust and accurate prediction of both binding poses and their affinities at early stage in designing drug candidates. Small, rigid ligands with few intramolecular degrees of freedom, for example, fragment-like molecules, have multiple binding poses, even at a single binding site, and their affinities are often close to each other. We explore various structures of ligand binding to a target through metadynamics using a small number of collective variables, followed by reweighting to obtain the atomic coordinates. After identifying each binding pose by cluster analysis, we perform alchemical free energy calculations on each structure to obtain the overall value. We applied this protocol in computing free energy of binding for the theophylline-RNA aptamer complex. Of the six (meta)stable structures found, the most favorable binding structure is consistent with the structure obtained by NMR. The overall free energy of binding reproduces the experimental values very well.

Modern drug design: the implication of using artificial neuronal networks and multiple molecular dynamic simulations

We report the implementation of molecular modeling approaches developed as a part of the 2016 Grand Challenge 2, the blinded competition of computer aided drug design technologies held by the D3R Drug Design Data Resource ( https://drugdesigndata.org/ ). The challenge was focused on the ligands of the farnesoid X receptor (FXR), a highly flexible nuclear receptor of the cholesterol derivative chenodeoxycholic acid. FXR is considered an important therapeutic target for metabolic, inflammatory, bowel and obesity related diseases (Expert Opin Drug Metab Toxicol 4:523-532, 2015), but in the context of this competition it is also interesting due to the significant ligand-induced conformational changes displayed by the protein. To deal with these conformational changes we employed multiple simulations of molecular dynamics (MD). Our MD-based protocols were top-ranked in estimating the free energy of binding of the ligands and FXR protein. Our approach was ranked second in the prediction of the binding poses where we also combined MD with molecular docking and artificial neural networks. Our approach showed mediocre results for high-throughput scoring of interactions.

Free Energy of Binding and Mechanism of Interaction for the MEEVD-TPR2A Peptide-Protein Complex

The association between the MEEVD C-terminal peptide from the heat shock protein 90 (Hsp90) and tetratricopeptide repeat A (TPR2A) domain of the heat shock organizing protein (Hop) is a useful prototype to study the fundamental molecular details about the Hop-Hsp90 interaction. We study here the mechanism of binding/unbinding and compute the standard binding free energy and potential of mean force for the association of the MEEVD peptide to the TPR2A domain using the Adaptive Biasing Force (ABF) methodology. We observe conformational changes of the peptide and the protein receptor induced by binding. We measure the binding free energy of -8.4 kcal/mol, which is consistent with experimental estimates. The simulations achieve multiple unbinding and rebinding events along a consistent pathway connecting the binding site to solvent. The MEEVD peptide slowly dissociates disrupting the hydrogen bonds first, then tilting on the side while preserving the interaction with the side chain of residue Asp 5 of the peptide. After this initial displacement, the peptide completely dissociates and moves into the solvent. Rebinding of the MEEVD peptide from the solvent to the receptor binding site occurs slowly through the portal of entry. Unbinding and rebinding go through intermediate states characterized by the peptide interacting with a lateral helix, helix A1, of the receptor with mainly Asp 5, Val 4, and Glu 3 of the peptide. This newly discovered intermediate structure is characterized by numerous contacts with the receptor which lead to complete formation of the bound complex. The structure of the bound complex obtained after rebinding is structurally very similar to the crystal structure of the complex (0.48 Å root-mean square deviation). The residues Asp 5, Val 4, and Glu 3 adopt conformations and intermolecular contacts with excellent structural similarity with the native ones. Finally, the dissociation and reassociation of MEEVD induce hydration/dehydration transitions, which provide insights on the role of desolvation and solvation processes in protein-peptide binding.

I. Dissociation free energies of drug–receptor systems via non-equilibrium alchemical simulations: a theoretical framework

In this contribution I critically discuss the alchemical approach for evaluating binding free energies in drug–receptor systems, placing this methodology into the broader context of non-equilibrium thermodynamics.

Evaluating the Accuracy of the Quasiharmonic Approximation

The quasiharmonic approximation (QH) allows the configurational entropy of a molecule to be estimated on the basis of a molecular dynamics simulation, through construction of a Gaussian probability distribution of conformations with variances equal to those provided by the simulation. At its introduction in 1981, the QH method was successfully applied to simple molecular systems with only one highly occupied energy well, and fluctuations were analyzed in a system of internal bond-angle-torsion coordinates. However, more recent studies have applied the QH method to complex biomolecular systems and have relied upon Cartesian coordinates. The present study evaluates the accuracy of the QH method through comparisons with more detailed methods. The chief findings are that the QH method can markedly overestimate the configurational entropy for systems with multiple occupied energy wells and that such errors tend to be magnified by the use of Cartesian coordinates instead of bond-angle-torsion coordinates.

Calculation of Standard Binding Free Energies: Aromatic Molecules in the T4 Lysozyme L99A Mutant

Calculations of the binding free energy of various nonpolar aromatic ligands with the L99A mutant of T4 lysozyme using molecular dynamics (MD) simulation are presented. To ensure better convergence, biasing potentials are used to restrain the ligand orientation and center-of-mass movement relative to the binding site when the ligand is decoupled from its environment in the binding pocket. The bias introduced by the restraint potentials is removed once the ligand fully interacts with the rest of the system and the calculated binding free energy is independent of the applied restraints. To decrease the computational cost, the simulations are generated with a reduced system in which protein and water atoms within a 15 Å-radius sphere around the ligand are included explicitly, while the rest of the system is treated with the generalized solvent boundary potential (GSBP). For all the ligands, the precision of the calculated free energy is less than 0.5 kcal/mol. For small nonpolar ligands such as benzene, toluene, and ethylbenzene, the calculated binding free energies are within 1.1 kcal/mol of the experimental values. For larger ligands, the computed binding free energies are slightly more favorable than the experimental values. The nonbinding polar molecule, phenol, has a calculated binding free energy of -0.88 kcal/mol. The simulation protocol presented here provides a way to calculate the binding free energy of small molecules to receptors at moderate computational cost.

Practical Aspects of Free-Energy Calculations: A Review

Free-energy calculations in the framework of classical molecular dynamics simulations are nowadays used in a wide range of research areas including solvation thermodynamics, molecular recognition, and protein folding. The basic components of a free-energy calculation, that is, a suitable model Hamiltonian, a sampling protocol, and an estimator for the free energy, are independent of the specific application. However, the attention that one has to pay to these components depends considerably on the specific application. Here, we review six different areas of application and discuss the relative importance of the three main components to provide the reader with an organigram and to make nonexperts aware of the many pitfalls present in free energy calculations.

Protein-ligand binding free energies from exhaustive docking

We explore the use of exhaustive docking as an alternative to Monte Carlo and molecular dynamics sampling for the direct integration of the partition function for protein-ligand binding. We enumerate feasible poses for the ligand and calculate the Boltzmann factor contribution of each pose to the partition function. From the partition function, the free energy, enthalpy, and entropy can be derived. All our calculations are done with a continuum solvation model that includes solving the Poisson equation. In contrast to Monte Carlo and molecular dynamics simulations, exhaustive docking avoids (within the limitations of a discrete sampling) the question of "Have we run long enough?" due to its deterministic complete enumeration of states. We tested the method on the T4 lysozyme L99A mutant, which has a nonpolar cavity that can accommodate a number of small molecules. We tested two electrostatic models. Model 1 used a solute dielectric of 2.25 for the complex apoprotein and free ligand and 78.5 for the solvent. Model 2 used a solute dielectric of 2.25 for the complex and apoprotein but 1.0 for the free ligand. For our test set of eight molecules, we obtain a reasonable correlation with a Pearson r(2) = 0.66 using model 1. The trend in binding affinity ranking is generally preserved with a Kendall τ = 0.64 and Spearman ρ = 0.83. With model 2, the correlation is improved with a Pearson r(2) = 0.83, Kendall τ = 0.93, and Spearman ρ = 0.98. This suggests that the energy function and sampling method adequately captured most of the thermodynamics of binding of the nonpolar ligands to T4 lysozyme L99A.

Blind prediction of host-guest binding affinities: a new SAMPL3 challenge

The computational prediction of protein-ligand binding affinities is of central interest in early-stage drug-discovery, and there is a widely recognized need for improved methods. Low molecular weight receptors and their ligands--i.e., host-guest systems--represent valuable test-beds for such affinity prediction methods, because their small size makes for fast calculations and relatively facile numerical convergence. The SAMPL3 community exercise included the first ever blind prediction challenge for host-guest binding affinities, through the incorporation of 11 new host-guest complexes. Ten participating research groups addressed this challenge with a variety of approaches. Statistical assessment indicates that, although most methods performed well at predicting some general trends in binding affinity, overall accuracy was not high, as all the methods suffered from either poor correlation or high RMS errors or both. There was no clear advantage in using explicit versus implicit solvent models, any particular force field, or any particular approach to conformational sampling. In a few cases, predictions using very similar energy models but different sampling and/or free-energy methods resulted in significantly different results. The protonation states of one host and some guest molecules emerged as key uncertainties beyond the choice of computational approach. The present results have implications for methods development and future blind prediction exercises.

Role of Ligand Reorganization and Conformational Restraints on the Binding Free Energies of DAPY Non-Nucleoside Inhibitors to HIV Reverse Transcriptase

The results of computer simulations of the binding of etravirine (TMC125) and rilpivirine (TMC278) to HIV reverse transcriptase are reported. It is confirmed that consistent binding free energy estimates are obtained with or without the application of torsional restraints when the free energies of imposing the restraints are taken into account. The restraints have a smaller influence on the thermodynamics and apparent kinetics of binding of TMC125 compared to the more flexible TMC278 inhibitor. The concept of the reorganization free energy of binding is useful to understand and categorize these effects. Contrary to expectations, the use of conformational restraints did not consistently enhance convergence of binding free energy estimates due to suppression of binding/unbinding pathways and due to the influence of rotational degrees of freedom not directly controlled by the restraints. Physical insights concerning the thermodynamic driving forces for binding and the role of "jiggling" and "wiggling" motion of the ligands are discussed. Based on these insights we conclude that an ideal inhibitor, if chemically realizable, would possess the electrostatic charge distribution of TMC125, so as to form strong interactions with the receptor, and the larger and more flexible substituents of TMC278, so as to minimize reorganization free energy penalties and the effects of resistance mutations, suitably modified, as in TMC125, so as to disfavor the formation of non-binding competent extended conformations when free in solution.

Conformational Transitions and Convergence of Absolute Binding Free Energy Calculations

The Binding Energy Distribution Analysis Method (BEDAM) is employed to compute the standard binding free energies of a series of ligands to a FK506 binding protein (FKBP12) with implicit solvation. Binding free energy estimates are in reasonably good agreement with experimental affinities. The conformations of the complexes identified by the simulations are in good agreement with crystallographic data, which was not used to restrain ligand orientations. The BEDAM method is based on λ -hopping Hamiltonian parallel Replica Exchange (HREM) molecular dynamics conformational sampling, the OPLS-AA/AGBNP2 effective potential, and multi-state free energy estimators (MBAR). Achieving converged and accurate results depends on all of these elements of the calculation. Convergence of the binding free energy is tied to the level of convergence of binding energy distributions at critical intermediate states where bound and unbound states are at equilibrium, and where the rate of binding/unbinding conformational transitions is maximal. This finding mirrors similar observations in the context of order/disorder transitions as for example in protein folding. Insights concerning the physical mechanism of ligand binding and unbinding are obtained. Convergence for the largest FK506 ligand is achieved only after imposing strict conformational restraints, which however require accurate prior structural knowledge of the structure of the complex. The analytical AGBNP2 model is found to underestimate the magnitude of the hydrophobic driving force towards binding in these systems characterized by loosely packed protein-ligand binding interfaces. Rescoring of the binding energies using a numerical surface area model corrects this deficiency. This study illustrates the complex interplay between energy models, exploration of conformational space, and free energy estimators needed to obtain robust estimates from binding free energy calculations.

Advances in all atom sampling methods for modeling protein-ligand binding affinities

Conformational dynamics plays a fundamental role in the regulation of molecular recognition processes. Conformational heterogeneity and entropy variations upon binding, although not always evident from the analysis of structural data, can substantially affect affinity and specificity. Computer modeling is able to provide some of the most direct insights into these aspects of molecular recognition. We review recent physics-based computational studies that employ advanced conformational sampling algorithms and effective potentials to model the three main classes of degrees of freedom relevant to the binding process: ligand positioning relative to the receptor, ligand and receptor internal reorganization, and hydration. Collectively these studies show that all of these elements are important for proper modeling of protein-ligand interactions.

Recent theoretical and computational advances for modeling protein-ligand binding affinities

We review recent theoretical and algorithmic advances for the modeling of protein ligand binding free energies. We first describe a statistical mechanics theory of noncovalent association, with particular focus on deriving the fundamental formulas on which computational methods are based. The second part reviews the main computational models and algorithms in current use or development, pointing out the relations with each other and with the theory developed in the first part. Particular emphasis is given to the modeling of conformational reorganization and entropic effect. The methods reviewed are free energy perturbation, double decoupling, the Binding Energy Distribution Analysis Method, the potential of mean force method, mining minima and MM/PBSA. These models have different features and limitations, and their ranges of applicability vary correspondingly. Yet their origins can all be traced back to a single fundamental theory.

Lipid-modulated sequence-specific association of glycophorin A in membranes

Protein association in lipid membranes is a complex process with thermodynamics directed by a multitude of different factors. Amino-acid sequence is a molecular parameter that affects dimerization as shown by limited directed mutations along the transmembrane domains. Membrane-mediated interactions are also important although details of such contributions remain largely unclear. In this study, we probe directly the free energy of association of Glycophorin A by means of extensive parallel Monte Carlo simulations with recently developed methods and a model that accounts for sequence-specificity while representing lipid membranes faithfully. We find that lipid-induced interactions are significant both at short and intermediate separations. The ability of molecules to tilt in a specific hydrophobic environment extends their accessible interfaces, leading to intermittent contacts during protein recognition. The dimer with the lowest free energy is largely determined by the favorable lipid-induced attractive interactions at the closest distance. Finally, the coarse-grained model employed herein, together with the extensive sampling performed, provides estimates of the free energy of association that are in excellent agreement with existing data.

The Binding Energy Distribution Analysis Method (BEDAM) for the Estimation of Protein-Ligand Binding Affinities

The Binding Energy Distribution Analysis Method (BEDAM) for the computation of receptor-ligand standard binding free energies with implicit solvation is presented. The method is based on a well established statistical mechanics theory of molecular association. It is shown that, in the context of implicit solvation, the theory is homologous to the test particle method of solvation thermodynamics with the solute-solvent potential represented by the effective binding energy of the protein-ligand complex. Accordingly, in BEDAM the binding constant is computed by means of a weighted integral of the probability distribution of the binding energy obtained in the canonical ensemble in which the ligand is positioned in the binding site but the receptor and the ligand interact only with the solvent continuum. It is shown that the binding energy distribution encodes all of the physical effects of binding. The balance between binding enthalpy and entropy is seen in our formalism as a balance between favorable and unfavorable binding modes which are coupled through the normalization of the binding energy distribution function. An efficient computational protocol for the binding energy distribution based on the AGBNP2 implicit solvent model, parallel Hamiltonian replica exchange sampling and histogram reweighting is developed. Applications of the method to a set of known binders and non-binders of the L99A and L99A/M102Q mutants of T4 lysozyme receptor are illustrated. The method is able to discriminate without error binders from non-binders, and the computed standard binding free energies of the binders are found to be in good agreement with experimental measurements. Analysis of the results reveals that the binding affinities of these systems reflect the contributions from multiple conformations spanning a wide range of binding energies.

A schematic model for molecular affinity and binding with Ising variables

After discussing the relevance of statistical physics in molecular recognition processes, we present a schematic model for ligand-receptor association based on an Ising chain. We discuss the possible behaviors of the affinity when the stiffness of the ligand increases. We also consider the case of flexible receptors. A variety of interesting behaviors is obtained, including some affinity modulation upon bond hardening or softening. The affinity of a ligand for its receptor is shown to depend on the details of its rigidity profile, and we question the possibility of encoding information in the rigidities as well as in the shape. An exhaustive study of the selectivity of patterns with length n < 8 is carried out. Connection with other spin models, in particular spin glasses is mentioned in the conclusion.

Entropic cost of protein-ligand binding and its dependence on the entropy in solution

Two theoretical formulations are proposed and compared for the loss of translational and rotational entropy upon protein-ligand binding in water. The two theories share the same approach to evaluate the translational and rotational entropy of the ligand when bound. The potential of the bound ligand is modeled by six harmonic oscillators that are parametrized from the force and torque magnitudes measured in a molecular dynamics simulation, yielding vibrational and librational entropies. In the aqueous phase, the theories differ because there is no unique way to assign the total entropy to molecules in solution. In one approach, the ligand is allowed unrestricted access to the full solution volume at the standard concentration and is assigned the same translational and rotational entropy as if it were an ideal gas. We term this a "molecule-frame" (MF) theory because it considers configurational space in the reference frame of the molecule of interest. The entropy of the solvent is penalized because it is excluded from the molecule's volume. In the second theory, all molecules including the solvent are confined by their neighbors in mean-field configurational volumes. This we term a "system-frame" (SF) theory because the configurational space available to all molecules is considered in the reference frame of the whole system. Molecules have vibrational and librational entropy in the same way as they do when bound. In addition, the discrete size of the solvent molecules quantizes the configurational space into an effective number of minima according to the solute molecule's standard concentration and the mean volume of a solvent molecule. This leads to the cratic entropy expressed in terms of the solute molecule's mole fraction. The equivalent number of minima in rotational space depends on both the solute molecule's volume and the solvent molecule's volume. This leads to an equation for the orientational entropy based on the proposed concept of "angle fraction". The MF and SF theories are applied to calculate the translational and rotational entropy losses involved in the formation of six different protein-ligand complexes, in two of which the ligand is water. The MF entropy losses range from -80 to -142 J K(-1) mol(-1) for ligands at the 1 M standard-state concentration and from -52 to -63 J K(-1) mol(-1) for water at the 55.6 M standard-state concentration. They depend logarithmically on both the number and strength of interactions between the ligand and protein through the forces and torques. This is observed to lead to moderate dependencies on the ligands' moments of inertia and masses. The SF entropy losses are smaller and range from -50 to -75 J K(-1) mol(-1) for ligands at the 1 M standard-state concentration and from 0 to -12 J K(-1) mol(-1) for water. They depend logarithmically on the ligand solvent's molecular volume and weakly on the relative strengths of the ligand's interactions with the protein and water. The cratic entropy loss in water at the standard concentration is constant and is also demonstrated to be implicit in MF theories. Entropy losses from the two approaches are also compared with those from other computational approaches and with experiment. The use of the force and torque magnitudes leads to smaller bound volumes than are obtained from ligand-displacement approaches. The general agreement of the SF entropy losses with those from experiment suggests that the SF theory is more consistent with the assumptions made in experimental measurements than the MF solvation theories, which would require a compensating entropy gain in the solvent in order to agree.

Standard Free Energy of Binding from a One-Dimensional Potential of Mean Force

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.

Computations of standard binding free energies with molecular dynamics simulations

An increasing number of studies have reported computations of the standard (absolute) binding free energy of small ligands to proteins using molecular dynamics (MD) simulations and explicit solvent molecules that are in good agreement with experiments. This encouraging progress suggests that physics-based approaches hold the promise of making important contributions to the process of drug discovery and optimization in the near future. Two types of approaches are principally used to compute binding free energies with MD simulations. The most widely known is the alchemical double decoupling method, in which the interaction of the ligand with its surroundings are progressively switched off. It is also possible to use a potential of mean force (PMF) method, in which the ligand is physically separated from the protein receptor. For both of these computational approaches, restraining potentials may be activated and released during the simulation for sampling efficiently the changes in translational, rotational, and conformational freedom of the ligand and protein upon binding. Because such restraining potentials add bias to the simulations, it is important that their effects be rigorously removed to yield a binding free energy that is properly unbiased with respect to the standard state. A review of recent results is presented, and differences in computational methods are discussed. Examples of computations with T4-lysozyme mutants, FKBP12, SH2 domain, and cytochrome P450 are discussed and compared. Remaining difficulties and challenges are highlighted.

Calculation of absolute ligand binding free energy to a ribosome-targeting protein as a function of solvent model

A comparative analysis is provided of the effect of different solvent models on the calculation of a potential of mean force (PMF) for determining the absolute binding affinity of the small molecule inhibitor pteroic acid bound to ricin toxin A-chain (RTA). Solvent models include the distance-dependent dielectric constant, several different generalized Born (GB) approximations, and a hybrid explicit/GB-based implicit solvent model. We found that the simpler approximation of dielectric screening and a GB model, with Born radii fitted to a switching-window dielectric-boundary surface Poisson solvent model, severely overpredicted the binding affinity as compared to the experimental value, estimated to range from -4.4 to -6.0 kcal/mol. In contrast, GB models that are parametrized to fit the Lee-Richards molecular surface performed much better, predicting binding free energy within 1-3 kcal/mol of experimental estimates. However, the predicted free-energy profiles of these GB models displayed alternative binding modes not observed in the crystal structure. Finally, the most rigorous and computationally costly approach in this work, which used a hybrid explicit/implicit solvent model, correctly determined a binding funnel in the PMF near the crystallographic bound state and predicted an absolute binding affinity that was 2 kcal/mol more favorable than the estimated experimental binding affinity.

Why are proteins charged? Networks of charge-charge interactions in proteins measured by charge ladders and capillary electrophoresis

Almost all proteins contain charged amino acids. While the function in catalysis or binding of individual charges in the active site can often be identified, it is less clear how to assign function to charges beyond this region. Are they necessary for solubility? For reasons other than solubility? Can manipulating these charges change the properties of proteins? A combination of capillary electrophoresis (CE) and protein charge ladders makes it possible to study the roles of charged residues on the surface of proteins outside the active site. This method involves chemical modification of those residues to generate a large number of derivatives of the protein that differ in charge. CE separates those derivatives into groups with the same number of modified charged groups. By studying the influence of charge on the properties of proteins using charge ladders, it is possible to estimate the net charge and hydrodynamic radius and to infer the role of charged residues in ligand binding and protein folding.

Calculation of Protein-Ligand Binding Affinities

Accurate methods of computing the affinity of a small molecule with a protein are needed to speed the discovery of new medications and biological probes. This paper reviews physics-based models of binding, beginning with a summary of the changes in potential energy, solvation energy, and configurational entropy that influence affinity, and a theoretical overview to frame the discussion of specific computational approaches. Important advances are reported in modeling protein-ligand energetics, such as the incorporation of electronic polarization and the use of quantum mechanical methods. Recent calculations suggest that changes in configurational entropy strongly oppose binding and must be included if accurate affinities are to be obtained. The linear interaction energy (LIE) and molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) methods are analyzed, as are free energy pathway methods, which show promise and may be ready for more extensive testing. Ultimately, major improvements in modeling accuracy will likely require advances on multiple fronts, as well as continued validation against experiment.

Binding pathways of ligands to HIV-1 protease: coarse-grained and atomistic simulations

Multiscale simulations (coarse-grained Brownian dynamics simulations and all-atom molecular dynamics simulations in implicit solvent) were applied to reveal the binding processes of ligands as they enter the binding site of the HIV-1 protease. The initial structures used for the molecular dynamics simulations were generated based on the Brownian dynamics trajectories, and this is the first molecular dynamics simulation of modeling the association of a ligand with the protease. We found that a protease substrate successfully binds to the protein when the flaps are fully open. Surprisingly, a smaller cyclic urea inhibitor (XK263) can reach the binding site when the flaps are not fully open. However, if the flaps are nearly closed, the inhibitor must rearrange or binding can fail because the inhibitor cannot attain proper conformations to enter the binding site. Both the peptide substrate and XK263 can also affect the protein's internal motion, which may help the flaps to open. Simulations allow us to efficiently study the ligand binding processes and may help those who study drug discovery to find optimal association pathways and to design those ligands with the best binding kinetics.

Energy landscape and transition state of protein-protein association

Formation of a stereospecific protein complex is favored by specific interactions between two proteins but disfavored by the loss of translational and rotational freedom. Echoing the protein folding process, we have previously proposed a transition state for protein-protein association. Here we clarify the specification of the transition state by working with two types of toy models for protein association. A "hemisphere" model consists of two matching hemispheres as associating proteins, and a "crater" model consists of a spherical protein with a crater to which another spherical protein fits snugly. Short-range pairwise interactions between sites across the interface hold together the bound complex. Small relative translation and rotation between the subunits quickly destroy the interactions, leading to a sharp transition between the bound state with numerous short-range interactions but restricted translation and rotational freedom and the unbound state with, at most, a small number of interactions but expanded configurational freedom. This transition sets the outer boundary of the bound state as well as the transition state for association. The energy landscape is funnel-like, with the deep well of the bound state surrounded by a broad shallow basin. Calculations with the toy models suggest that mutational change in the interaction energy in the x-ray structure of a protein-protein complex, commonly used to approximate the effect on the association constant, overestimates the effect of mutation by 10-20%. With an eye toward specifying the transition states of actual protein complexes, we find that the total number of contacts between the subunits serves as a good surrogate of the interaction energy. This formalism of protein-protein association is applied to the barnase-barstar complex, reproducing the experimental results for the association rate over a wide range of ionic strength.

Analysis of the conformational dependence of mass-metric tensor determinants in serial polymers with constraints

It is well known that mass-metric tensor determinants det(G(s)) influence the equilibrium statistics and the rates of conformational transitions for polymers with constrained bond lengths and bond angles. It is now standard practice to include a Fixman-style compensating potential of the form U(c)(q(s)) proportional, variant(-k(B)T/2)ln[det(G(s))] as part of algorithms for torsional space molecular dynamics. This elegant strategy helps eliminate unwarranted biases that arise due to the imposition of holonomic constraints. However, the precise nature and extent of variation of det(G(s)) and hence ln[det(G(s))] with chain conformation and chain length has never been quantified. This type of analysis is crucial for understanding the nature of the conformational bias that the introduction of a Fixman potential aims to eliminate. Additionally, a detailed analysis of the conformational dependence of det(G(s)) will help resolve ambiguities regarding suggestions for incorporating terms related to det(G(s)) in the design of move sets in torsional space Monte Carlo simulations. In this work, we present results from a systematic study of the variation of det(G(s)) for a serial polymer with fixed bond lengths and bond angles as a function of chain conformation and chain length. This analysis requires an algorithm designed for rapid computation of det(G(s)) which simultaneously allows for a physical/geometric interpretation of the conformational dependence of det(G(s)). Consequently, we provide a detailed discussion of our adaptation of an O(n) algorithm from the robotics literature, which leads to simple recursion relations for direct evaluation of det(G(s)). Our analysis of the conformational dependence of det(G(s)) yields the following insights. (1) det(G(s)) is maximized for spatial conformers and minimized for planar conformations. (2) Previous work suggests that it is logical to expect that the conformational dependence of det(G(s)) becomes more pronounced with increase in chain length. Confirming this expectation, we provide systematic quantification of the nature of this dependency and show that the difference in det(G(s)) between spatial and planar conformers, i.e., between the maxima and minima of det(G(s)) grows systematically with chain length. Finally, we provide a brief discussion of implications of our analysis for the design of move sets in Monte Carlo simulations.

Calculation of absolute protein-ligand binding affinity using path and endpoint approaches

A comparative analysis is provided of rigorous and approximate methods for calculating absolute binding affinities of two protein-ligand complexes: the FKBP protein bound with small molecules 4-hydroxy-2-butanone and FK506. Our rigorous approach is an umbrella sampling technique where a potential of mean force is determined by pulling the ligand out of the protein active site over several simulation windows. The results of this approach agree well with experimentally observed binding affinities. Also assessed is a commonly used approximate endpoint approach, which separately estimates enthalpy, solvation free energy, and entropy. We show that this endpoint approach has numerous variations, all of which are prone to critical shortcomings. For example, conventional harmonic and quasiharmonic entropy estimation procedures produce disparate results for the relatively simple protein-ligand systems studied in this work.

Free energy, entropy, and induced fit in host-guest recognition: calculations with the second-generation mining minima algorithm

This study applies a novel computational method to study molecular recognition for three sets of synthetic hosts: molecular clips, molecular tweezers, and a synthetic barbiturate receptor. The computed standard free energies of binding for the 12 binding reactions agree closely with experiment and provide insight into the roles of configurational entropy, preorganization, and induced fit in the systems studied. The computed changes in configurational entropy are comparable in magnitude to the changes in mean potential plus solvation energy, and they result primarily from changes in the average width of the energy wells upon binding. A strong correlation is observed between the changes in configurational energy and configurational entropy upon binding, resulting in near-linear compensation analogous to classical entropy-enthalpy compensation.