Absolute and relative binding free energy calculations of the interaction of biotin and its analogs with streptavidin using molecular dynamics/free energy perturbation approaches

The statistical-thermodynamic basis for computation of binding affinities: a critical review

Although the statistical thermodynamics of noncovalent binding has been considered in a number of theoretical papers, few methods of computing binding affinities are derived explicitly from this underlying theory. This has contributed to uncertainty and controversy in certain areas. This article therefore reviews and extends the connections of some important computational methods with the underlying statistical thermodynamics. A derivation of the standard free energy of binding forms the basis of this review. This derivation should be useful in formulating novel computational methods for predicting binding affinities. It also permits several important points to be established. For example, it is found that the double-annihilation method of computing binding energy does not yield the standard free energy of binding, but can be modified to yield this quantity. The derivation also makes it possible to define clearly the changes in translational, rotational, configurational, and solvent entropy upon binding. It is argued that molecular mass has a negligible effect upon the standard free energy of binding for biomolecular systems, and that the cratic entropy defined by Gurney is not a useful concept. In addition, the use of continuum models of the solvent in binding calculations is reviewed, and a formalism is presented for incorporating a limited number of solvent molecules explicitly.

Protein-biotin interactions

The three-dimensional structures of several biotin-binding proteins are now known, giving insights into the molecular architecture of the binding sites for biotin. In combination with biochemical and computational approaches, these structural insights provide the basis for our present understanding of biotin-protein interactions which, in some cases, give rise to spectacular binding constants.

On the bioactive conformation of a small peptide and its set of thermodynamically accessible conformations

In order to investigate the relationship between the bioactive conformation of a peptide and its set of thermodynamically accessible structures in solution, the conformational profile of the tetrapeptide Ac-Pro-Ala-Pro-Tyr-OH was characterized by computational methods. Search of the conformational space was performed within the molecular mechanics frame-work using the AMBER4.0 force field with an effective dielectric constant of 80. Unique structures of the peptide were compared with its bioactive conformation for the protein Streptomyces griseus Protease A, as taken from the crystal structure of the enzyme-peptide complex. The results show that the bound conformation is close to one of the unique conformations characterized in the conformational search of the isolated peptide. Moreover, the lowest energy minimum characterized in the conformational search exhibits large deviations when compared to the bound conformation of the crystal structure.

Scoring noncovalent protein-ligand interactions: a continuous differentiable function tuned to compute binding affinities

Exploitation of protein structures for potential drug leads by molecular docking is critically dependent on methods for scoring putative protein-ligand interactions. An ideal function for scoring must exhibit predictive accuracy and high computational speed, and must be tolerant of variations in the relative protein-ligand molecular alignment and conformation. This paper describes the development of an empirically derived scoring function, based on the binding affinities of protein-ligand complexes coupled with their crystallographically determined structures. The function's primary terms involve hydrophobic and polar complementarity, with additional terms for entropic and solvation effects. The issue of alignment/conformation dependence was solved by constructing a continuous differentiable nonlinear function with the requirement that maxima in ligand conformation/alignment space corresponded closely to crystallographically determined structures. The expected error in the predicted affinity based on cross-validation was 1.0 log unit. The function is sufficiently fast and accurate to serve as the objective function of a molecular-docking search engine. The function is particularly well suited to the docking problem, since it has spatially narrow maxima that are broadly accessible via gradient descent.

The enthalpy change in protein folding and binding: refinement of parameters for structure-based calculations

Two effects are mainly responsible for the observed enthalpy change in protein unfolding: the disruption of internal interactions within the protein molecule (van der Waals, hydrogen bonds, etc.) and the hydration of the groups that are buried in the native state and become exposed to the solvent on unfolding. In the traditional thermodynamic analysis, the effects of hydration have usually been evaluated using the thermodynamic data for the transfer of small model compounds from the gas phase to water. The contribution of internal interactions, on the other hand, are usually estimated by subtracting the hydration effects from the experimental enthalpy of unfolding. The main drawback of this approach is that the enthalpic contributions of hydration, and those due to the disruption of internal interactions, are more than one order of magnitude larger than the experimental enthalpy value. The enthalpy contributions of hydration and disruption of internal interactions have opposite signs and cancel each other almost completely resulting in a final value that is over 10 times smaller than the individual terms. For this reason, the classical approach cannot be used to accurately predict unfolding enthalpies from structure: any error in the estimation of the hydration enthalpy will be amplified by a factor of 10 or more in the estimation of the unfolding enthalpy. Recently, it has been shown that simple parametric equations that relate the enthalpy change with certain structural parameters, especially changes in solvent accessible surface areas have considerable predictive power. In this paper, we provide a physical foundation to that parametrization and in the process we present a system of equations that explicitly includes the enthalpic effects of the packing density between the different atoms within the protein molecule. Using this approach, the error in the prediction of folding/unfolding enthalpies at 60 degrees C, the median temperature for thermal unfolding, is better than +/- 3% (standard deviation = 4 kcal/mol).

Molecular docking using surface complementarity

A method is described to dock a ligand into a binding site in a protein on the basis of the complementarity of the intermolecular atomic contacts. Docking is performed by maximization of a complementarity function that is dependent on atomic contact surface area and the chemical properties of the contacting atoms. The generality and simplicity of the complementarity function ensure that a wide range of chemical structures can be handled. The ligand and the protein are treated as rigid bodies, but displacement of a small number of residues lining the ligand binding site can be taken into account. The method can assist in the design of improved ligands by indicating what changes in complementarity may occur as a result of the substitution of an atom in the ligand. The capabilities of the method are demonstrated by application to 14 protein-ligand complexes of known crystal structure.

Leucine 18, a hydrophobic residue essential for high affinity binding of anthopleurin B to the voltage-sensitive sodium channel

Anthopleurin B is a potent anemone toxin that binds with nanomolar affinity to the cardiac and neuronal isoforms of the voltage-gated sodium channel. A cationic cluster that includes Arg-12, Arg-14 and Lys-49 has been shown previously to be important in this interaction. In this study, we have used site-directed mutagenesis to determine the contribution to activity of two aliphatic residues, Leu-18 and Ile-43, that have previously been experimentally inaccessible. Leu-18, a residue proximal to the cationic cluster, plays a critical role in defining the high affinity of the toxin. In ion flux studies, this is exemplified by the several hundredfold loss in affinity (231-672-fold) observed for both L18A and L18V toxins on either isoform of the sodium channel. When analyzed electrophysiologically, L18A, the most severely compromised mutant, also displays a substantial loss in affinity (34-fold and 328-fold) for the neuronal and cardiac isoforms. This difference in affinities may reflect an increased preference of the L18A mutant for the closed state of the neuronal channel. In contrast, Ile-43, a residue distal to the cationic cluster, plays at most a very modest role in affinity toward both isoforms of the sodium channel. Only conservative substitutions are tolerated at this position, implying that it may contribute to an important structural component. Our results indicate that Leu-18 is the most significant single contributor to the high affinity of Anthopleurin B identified to date. These results have extended the binding site beyond the cationic cluster to include Leu-18 and broadened our emphasis from the basic residues to include the crucial role of hydrophobic residues in toxin-receptor interactions.

A technique to study molecular recognition in drug design: preliminary application of free energy derivatives to inhibition of a malarial cysteine protease

We present molecular dynamics studies on model complexes of inhibitors of a malarial cysteine protease. The initial model for such complexes came from the model building of the protein using its homology with other cysteine proteases and calculations using DOCK to generate new lead compounds. Some of the initial model-built structures were quite stable for 100 psec of dynamics; others moved significantly from their model-built orientation. We also calculated the free energy derivatives at each atom in the inhibitor, both in water and in the binding site. The results of these calculations suggest directions for the design of new, more potent enzyme inhibitors and agree qualitatively with some of the experimental findings. Nonetheless, we stress that we have only used this methodology in an interpretive rather than a predictive manner.

Ligand binding: molecular mechanics calculation of the streptavidin-biotin rupture force

The force required to rupture the streptavidin-biotin complex was calculated here by computer simulations. The computed force agrees well with that obtained by recent single molecule atomic force microscope experiments. These simulations suggest a detailed multiple-pathway rupture mechanism involving five major unbinding steps. Binding forces and specificity are attributed to a hydrogen bond network between the biotin ligand and residues within the binding pocket of streptavidin. During rupture, additional water bridges substantially enhance the stability of the complex and even dominate the binding interactions. In contrast, steric restraints do not appear to contribute to the binding forces, although conformational motions were observed.

Recombinant core streptavidins. A minimum-sized core streptavidin has enhanced structural stability and higher accessibility to biotinylated macromolecules

Two recombinant core streptavidins were designed and characterized to understand the role of the terminal sequences, present in naturally truncated core streptavidins, in the properties of streptavidin. One recombinant core streptavidin, Stv-25, has an amino acid sequence very similar to natural core streptavidins. The other recombinant molecule, Stv-13, has further truncation of the terminal residues and consists essentially of only the beta-barrel structure characteristic of streptavidin. These recombinant core streptavidins are tetrameric and bind four biotins/molecule, as does natural streptavidin. The solubility characteristics of Stv-13, determined by varying the concentration of ammonium sulfate or ethanol, were almost the same as those of Stv-25 and natural core streptavidin. However, Stv-13 showed an enhanced structural stability compared with Stv-25 and natural core streptavidin. For example, Stv-13 retained greater than 80% of its biotin binding ability after incubation in 6 M guanidine hydrochloride at pH 1.5, under which conditions, Stv-25 and natural core streptavidin retained only about 20% of their biotin binding ability. In addition, Stv-13 showed higher accessibility to biotinylated DNA than natural core streptavidin. Apparently, the terminal regions, present on the surface of natural core streptavidin, spatially hinder biotinylated macromolecules from approaching the biotin binding sites.

The relationship between ligand-binding thermodynamics and protein-ligand interaction forces measured by atomic force microscopy

The interaction forces between biotin and a set of streptavidin site-directed mutants with altered biotin-binding equilibrium and activation thermodynamics have been measured by atomic force microscopy. The AFM technique readily discriminates differences in interaction force between the site-directed (Trp to Phe or Ala) mutants. The interaction force is poorly correlated with both the equilibrium free energy of biotin binding and the activation free energy barrier to dissociation of the biotin-streptavidin complex. The interaction force is generally well correlated with the equilibrium biotin-binding enthalpy as well as the enthalpic activation barrier, but in the one mutant where these two parameters are altered in opposite directions, the interaction force is clearly correlated with the activation enthalpy of dissociation. These results suggest that the AFM force measurements directly probe the enthalpic activation barrier to ligand dissociation.

Intersubunit contacts made by tryptophan 120 with biotin are essential for both strong biotin binding and biotin-induced tighter subunit association of streptavidin

In natural streptavidin, tryptophan 120 of each subunit makes contacts with the biotin bound by an adjacent subunit through the dimer-dimer interface. To understand quantitatively the role of tryptophan 120 and its intersubunit communication in the properties of streptavidin, a streptavidin mutant in which tryptophan 120 is converted to phenylalanine was produced and characterized. The streptavidin mutant forms a tetrameric molecule and binds one biotin per subunit, as does natural streptavidin, indicating that the mutation of tryptophan 120 to phenylalanine has no significant effect on the basic properties of streptavidin. However, its biotin-binding affinity was reduced substantially, to approximately 10(8) M-1, indicating that the contact made by tryptophan 120 to biotin has a considerable contribution to the extremely tight biotin binding by streptavidin. The mutant retained bound biotin over a wide pH range or with the addition of urea up to 6 M at neutral pH. However, bound biotin was efficiently released by the addition of excess free biotin due, presumably, to exchange reactions. Electrophoretic analysis revealed that the intersubunit contact made by tryptophan 120 to biotin through the dimer-dimer interface is the major interaction responsible for the biotin-induced, tighter subunit association of streptavidin. In addition, the mutant has weaker subunit association than natural streptavidin even in the absence of biotin, indicating that tryptophan 120 also contributes to the subunit association of tetramers in the absence of biotin.

Protein motion and lock and key complementarity in antigen-antibody reactions

Antibodies possess a highly complementary combining site structure to that of their specific antigens. In many instances their reactions are driven by enthalpic factors including, at least in the case of the reaction of monoclonal antibody D1.3 with lysozyme, enthalpy of solvation. They require minor structural rearrangements, and their equilibrium association constants are relatively high (10(7)-10(11) M-1). By contrast, in an idiotope--anti-idiotope (antibody-antibody) reaction, which is entropically driven, the binding equilibrium constant is only 1.5 x 10(5) M-1 at 20 degrees C. This low value results from a slow association rate (10(3) M-1 s-1) due to a selection of conformational states that allow one of the interacting molecular surfaces (the idiotope on antibody D1.3) to become complementary to that of the anti-idiotopic antibody. Thus, antibody D1.3 reacts with two different macromolecules: with its specific antigen, hen egg lysozyme, and with a specific anti-idiotopic antibody. Complementarity with lysozyme is closer to a "lock and key" model and results in high affinity (2-4 x 10(8) M-1). That with the anti-idiotopic antibody involves conformational changes at its combining site and it results in a lower association constant (1.5 x 10(5) M-1). Thus, an "induced fit" mechanism may lead to a broadening of the binding specificity but with a resulting decrease in the intrinsic binding affinity which may weaken the physiological function of antibodies.

Site-directed mutagenesis studies of the high-affinity streptavidin-biotin complex: contributions of tryptophan residues 79, 108, and 120

We report the functional characterization of site-directed biotin binding-site mutants of recombinant core streptavidin. The mutagenesis studies were aimed at characterizing the contributions of Trp residues known to contact biotin that have been postulated to control the exceptional binding affinity observed in this system. The functional properties of single site-directed mutants replacing Trp residues with Phe or Ala at positions 79, 108, and 120 were investigated by quantitating the EC50 binding parameters of these mutants to biotin and 2-iminobiotin in an ELISA format. The biotin EC50 for all mutants was the same as wild-type streptavidin, demonstrating that their delta Ka values relative to wild type were < 10(6). The conservative W79F and W108F mutants displayed only a 2- to 3-fold increase in EC50 for 2-iminobiotin, corresponding to an estimated delta Ka < 10, while the W120F mutant displayed a much greater alteration in 2-iminobiotin EC50, corresponding to an estimated delta Ka of 10(2). These delta Ka values are likely to reflect similar changes for biotin. The 2-iminobiotin EC50 values for the Ala mutants fell outside the accessible concentration range of the ELISA assay, demonstrating that these mutations lowered the Ka by a factor of 10(4) to 10(6). Direct estimation of biotin Ka values for W79A, W120A, and W120F in an ultrafiltration binding assay yielded Ka values of 4.3 x 10(7) M-1, 8.6 x 10(6) M-1, and > 5 x 10(9) M-1, respectively, in excellent agreement with the ELISA estimates of delta Ka with 2-iminobiotin as a reporter ligand. The results of these preliminary functional studies suggest that these aromatic side chains contribute significantly to the streptavidin-biotin binding free energy.

Elusive affinities

The affinity of two molecules for each other and its temperature dependence are determined by the change in enthalpy, free enthalpy, entropy, and heat capacity upon dissociation. As we know the forces that stabilize protein-protein or protein-DNA association and the three-dimensional structures of the complex, we can in principle derive values for each one of these parameters. The calculation is done first in gas phase by molecular mechanics, then in solution with the help of hydration parameters calibrated on small molecules. However, estimates of enthalpy and entropy changes in gas phase have excessively large error bars even under the approximation that the components of the complex associate as rigid bodies. No reliable result can be expected at the end. The fit to experimental values derived from binding and calorimetric measurements is poor, except for the dissociation heat capacity. This parameter can be attributed mostly to the hydration step and it correlates with the size of the interface. Many protein-protein complexes have interface areas in the range 1200-2000 A2 and only small conformation changes, so the rigid body approximation applies. It is less generally valid in protein-DNA complexes, which have interfaces covering 2200-3100 A2, large dissociation heat capacities, and affect both the conformation and the dynamics of their components.

Model assembly study of the ligand binding by p-hydroxybenzoate hydroxylase: correlation between the calculated binding energies and the experimental dissociation constants

The energies of binding of seven ligands by p-hydroxybenzoate hydroxylase (PHBH) were calculated theoretically. Direct enzyme-ligand interaction energies were calculated using the ab initio quantum mechanical model assembly of the active site at the 3-21G level. Solvation energies of the ligands needed in the evaluation of the binding energies were calculated with the semiempirical AM1-SM2 method and the long-range electrostatic interaction energies between the ligands and the protein matrix classically using the static charge distributions of the ligands and the protein. Energies for proton-transfer between the ligands' OH or SH substituent at position 4 and the active-site tyrosine within the ab initio model assemblies were calculated and compared to the corresponding pKas in aqueous solution. Excluding 3,4-dihydroxybenzoate, the natural product of PHBH, a linear relationship between the calculated binding energies and the experimental binding free energies was found with a correlation coefficient of 0.90. Contributions of the direct enzyme-ligand interaction energies, solvation energies and the long-range electrostatic interaction energies to the calculated binding energies were analyzed. The proton-transfer energies of the ligands with substituents ortho to the ionized OH were found to be perturbed less in the model calculations than the energies of their meta isomers as deduced from the corresponding pKas.

Fundamentals of drug design from a biophysical viewpoint

Drug design means many things to many people. Commercially the aim is the development of compounds that can be patented and meet a variety of regulatory standards. In drug design, for medical purposes, toxicity and bio-availability are major considerations.

Quantum mechanical model assembly study on the energetics of binding of arabinose, fucose, and galactose to L-arabinose-binding protein

Binding energies of L-arabinose, D-fucose, and D-galactose to L-arabinose-binding protein was investigated theoretically. The calculated binding energies were composed of three contributions: 1) direct ligand-active site interaction energies calculated using static ab initio model assemblies; 2) solvation energies of the ligands; and 3) long-range electrostatic interaction energies between the ligands and the protein matrix. The calculated binding energies and the contributions of the energy components were used to analyze the experimental affinities of the ligands.

Intermolecular forces and energies between ligands and receptors

The recognition mechanisms and dissociation pathways of the avidin-biotin complex and of actin monomers in actin filaments were investigated. The unbinding forces of discrete complexes of avidin or streptavidin with biotin analogs are proportional to the enthalpy change of the complex formation but independent of changes in the free energy. This result indicates that the unbinding process is adiabatic and that entropic changes occur after unbinding. On the basis of the measured forces and binding energies, an effective rupture length of 9.5 +/- 1 angstroms was calculated for all biotin-avidin pairs and approximately 1 to 3 angstroms for the actin monomer-monomer interaction. A model for the correlation among binding forces, intermolecular potential, and molecular function is proposed.

Free energy simulations: the meaning of the individual contributions from a component analysis

A theoretical analysis is made of the decomposition into contributions from individual interactions of the free energy calculated by thermodynamic integration. It is demonstrated that such a decomposition, often referred to as "component analysis," is meaningful, even though it is a function of the integration path. Moreover, it is shown that the path dependence can be used to determine the relation of the contribution of a given interaction to the state of the system. To illustrate these conclusions, a simple transformation (Cl- to Br- in aqueous solution) is analyzed by use of the Reference Interaction Site Model-Hypernetted Chain Closure integral equation approach; it avoids the calculational difficulties of macromolecular simulation while retaining their conceptual complexity. The difference in the solvation free energy between chloride and bromide is calculated, and the contributions of the Lennard-Jones and electrostatic terms in the potential function are analyzed by the use of suitably chosen integration paths. The model is also used to examine the path dependence of individual contributions to the double free energy differences (delta delta G or delta delta A) that are often employed in free energy simulations of biological systems. The alchemical path, as contrasted with the experimental path, is shown to be appropriate for interpreting the effects of mutations on ligand binding and protein stability. The formulation is used to obtain a better understanding of the success of the Poisson-Boltzmann continuum approach for determining the solvation properties of polar and ionic systems.

Helix-capping interaction in lambda Cro protein: a free energy simulation analysis

The stability mutant Tyr-26-->Asp was studied in the Cro protein from bacteriophage lambda using free energy molecular dynamics simulations. The mutant was calculated to be more stable than the wild type by 3.0 +/- 1.7 kcal/mol/monomer, in reasonable agreement with experiment (1.4 kcal/mol/monomer). Moreover, the aspartic acid in the mutant was found to form a capping interaction with the amino terminus of the third alpha-helix of Cro. The simulations were analyzed to understand better the source of the stability of this helix-capping interaction and to examine the results in light of previous explanations of stabilizing helix caps--namely, a model of local unsatisfied hydrogen bonds at the helix termini and the helix macrodipole model. Analysis of the simulations shows that the stabilizing effect of this charged helical cap is due both to favorable hydrogen bonds with backbone NH groups at the helix terminus and to favorable electrostatic interactions (but not hydrogen bonds) with their carbonyls (effectively the next row of local dipoles in the helix). However, electrostatic interactions are weak or negligible with backbone dipolar groups in the helix further away from the terminus. Moreover, the importance of other local electrostatic interactions with polar side chains near the helix terminus, which are neglected in most treatments of this effect, are shown to be important. Thus, the results support a model that is intermediate between the two previous explanations: both unsatisfied hydrogen bonds at the helix terminus and other, local preoriented dipolar groups stabilize the helix cap. These findings suggest that similar interactions with preoriented dipolar groups may be important for cooperativity in other charge-dipole interactions and may be employed to advantage for molecular design.

What determines the strength of noncovalent association of ligands to proteins in aqueous solution?

Free energy perturbation methods using molecular dynamics have been used to calculate the absolute free energy of association of two ligand-protein complexes. The calculations reproduce the significantly more negative free energy of association of biotin to streptavidin, compared to N-L-acetyltryptophanamide/alpha-chymotrypsin. This difference in free energy of association is due to van der Waals/dispersion effects in the nearly ideally performed cavity that streptavidin presents to biotin, which involves four tryptophan residues.