INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

To understand molecular recognition of protein surfaces, it is important to obtain the quantitative contributions of the individual forces governing binding affinity and specificity. Consequently, it is useful to develop reliable methods for computing absolute or relative free energies that can dissect the relative magnitudes of different thermodynamics effects. Many such methods have been developed and they fall generally into three classes. The first class includes free energy simulation methods (Kollman, 1996; Tembe and McCammon, 1984) that treat solute and solvent atoms explicitly with an atomic force-field description. These simulation methods have a rigorous statistical mechanics basis and can be used to compute the difference between the binding free energies of two closely related medium-sized compounds (e.g., ligands differing only in a single residue) in solution fairly accurately (Jorgensen and Tirado-Rives, 1995). However, they are less reliable for ligand-protein or protein-protein association due mainly to the problem of sampling the relevant configurations of the solute and solvent (McCammon, 1998).

The second class encompasses empirical approaches that use parameters obtained from a "training set" of known interactions (Andrews et al., 1984; Bohm, 1994; Horton and Lewis, 1992; Weng et al., 1997; Williams et al., 1993). The binding free energy is taken to be a sum of different terms associated with hydrophobic contacts, hydrogen bonds or salt bridges, and conformational entropy loss. The various terms are usually not derived from theoretical backgrounds (e.g., statistical thermodynamics or molecular mechanics), but from available binding data in a statistical manner. The key advantage of empirical methods is their speed. However, these phenomenological scoring functions rely critically on the quality of the training database.

The third class includes methods that treat the solute atoms explicitly and the solvent as a continuum dielectric medium (Kollman et al., 2000). The electrostatic contribution to the solvation free energy is computed using finite difference Poisson-Boltzmann (PB) (Honig and Nicholls, 1995) methods or the Generalized Born (GB) approximation (Still et al., 1990; Qiu et al., 1997), whereas the nonelectrostatic contribution is treated as an empirical function of the solvent-accessible surface area (SASA) (Still et al., 1990). These methods appear promising for certain systems because they represent a good balance between speed (by avoiding sampling of solvent configurations) and accuracy (by incorporating long-range electrostatics, ionic strength, and polarization effects).

Although the computed absolute/relative binding free energies for various systems have been reported to be in close agreement with the respective experimental values, this does not necessarily mean that the underlying theory or scheme for computing the binding free energy is accurate. This is because the agreement may be fortuitous and errors due to experimental issues may mask the errors inherent in computing the various free energy components. For example, in many cases the x-ray structure of the complex, but not its unbound counterparts, is known. In such cases, the binding free energy is often estimated from the atomic coordinates of the complex alone assuming no conformational change upon complexation (rigid-binding approximation) (Novotny et al., 1989; Williams et al., 1991; Horton and Lewis, 1992; Jackson and Sternberg, 1995; Friedman and Honig, 1995; Froloff et al., 1997; Novotny et al., 1997; Olson, 1999). In cases where the atomic coordinates of both free and bound proteins are available, the protein may have been crystallized at a pH that is significantly different from the pH at which the binding free energy was measured. Furthermore, NMR structures are often solved at low pH. However, the free energy for protein association that is mediated by ionizable residues can be very sensitive to pH (Gibas, 1997; Xavier and Wilson, 1998). Discrepancy between the computed and experimental binding free energies has been attributed to the use of low pH structures rather than shortcomings in the free energy function/calculations, like the neglect of van der Waals (vdW) interactions (Krystek et al., 1993).

Partly due to the above issues, the various theoretical studies have not reached a consensus regarding the major force driving the binding process. Pauling concluded that the cooperation of weak vdW and hydrogen-bonding interactions between complementary surfaces over an area sufficiently large to overcome the disrupting influence of thermal agitation govern protein-protein recognition (Pauling, 1974; Pauling and Pressman, 1945). In contrast, Chothia and Janin suggested that vdW and polar interactions contribute little to complex stability, and it is hydrophobicity (the release of interfacial water) that is the major factor stabilizing protein-protein associations (Chothia and Janin, 1975; Kauzmann, 1959). This view is also implicitly assumed in studies that fit the experimental binding free energy to buried surface areas (Horton and Lewis, 1992). Honig and coworkers and other groups found that nonpolar interactions, represented by a free energy-surface area relationship that is assumed to account for the hydrophobic effect and enhanced interfacial packing, provide the major driving force for MHC class I protein-peptide, trypsin-inhibitor and cl repressor protein-DNA association (Froloff et al., 1997; Jayaram et al., 1999; Misra et al., 1998). Note that Chothia and Janin as well as Honig and coworkers did not take vdW interactions into account explicitly in computing the absolute binding free energy.

In this work, Scheme 1 (see next section), in conjunction with continuum dielectric models, has been implemented to compute the absolute free energy of protein-protein association. Our initial goals are twofold: to assess the accuracy of our scheme for computing absolute binding free energies, and to assess the validity of the various assumptions in previous theoretical treatments like the rigid-binding approximation or the neglect of vdW interactions and vibrational entropy (Novotny et al., 1989; Williams et al., 1991; Horton and Lewis, 1992; Vajda et al., 1994; Jackson and Sternberg, 1995; Novotny et al., 1997; Froloff et al., 1997; Olson, 1999). To address our first goal, we chose systems for which experimental binding free energies are available for comparison with the computed values, and accurate structures are available to minimize errors due to experimental issues. In this way, any observed discrepancy between the calculated and experimental results would reflect primarily the assumptions or deficiencies in our methodology rather than inaccuracies in the experimental structures. Specifically, we chose systems for which the crystal structures of both the complex and its components have been solved to 1.9-Å resolution to eliminate errors arising from poor quality structures or absence of free protein structures. Having the x-ray structures of the unbound and bound proteins also enable us to account for conformational changes that accompany binding (especially at the interaction surfaces). Also, we verified that the pH at which the free and complex structures were solved is close to the pH at which the corresponding binding free energy was measured to minimize errors due to structural changes at different pH values. To address our second goal, the systems chosen satisfied not only the above criteria, but have also been studied using various approximations.

Only two systems were found to satisfy all of the aforementioned criteria. They are trypsin complexed with bovine pancreatic trypsin inhibitor (BPTI), and the fragment variable (Fv) region of mouse monoclonal antibody, D1.3, bound to hen egg-white lysozyme (HEL). Their experimental binding free energies, Protein Data Bank (PDB) entries (Bernstein et al., 1977) and various properties are listed in Table 1. The interface area is defined as the SASA difference between the complex and its components. Note that both systems correspond to association of positively charged proteins to form a +12e complex with only small changes in conformation, as evidenced by root-mean-square deviation (RMSD) of the backbone atoms in the free x-ray structure from that in the corresponding complex crystal structure of <1 Š(see Table 1). The excellent agreement found between the computed and experimental binding free energies (see Results) allowed us to identify the key forces governing protease-inhibitor and antibody-antigen complexation with interface areas in the range of 1400-1700 Ų. Finally, we compare our findings with results obtained in previous studies for the same or homologous systems as well as for other protein-protein complexes (see Discussion).

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Theory

The binding free energy calculations were based on the thermodynamic cycle
In Scheme 1, lower case r and l denote the unbound receptor and ligand conformations, whereas upper case R and L denote the bound receptor and ligand conformations. Note that [x]_sln and [x]_gas, where x = r, l, or R · L, correspond to the same structure; viz., the relaxed "all-hydrogen" x-ray structure (see next section). The latter was energy minimized in the gas phase to yield [x^min]_gas, which corresponds to a local energy minimum that is required to compute gas-phase vibrational frequencies and thus the vibrational entropy change accompanying binding (see Gas-phase Binding Free Energy).

The standard free energy change (G⁰) for the noncovalent association of two molecules in solution according to Scheme 1 is given by

(1)

where

(2)

and

(3)

The desolvation free energy, G_solv, corresponds to the work of transferring the molecule in its solution conformation to the same conformation in the gas phase at 298 K. The isomerization free energy, G_isom, is the work of transforming the molecule in its x-ray conformation to a nearby local minimum structure in the gas phase. G_gas is the standard free energy change per mole for the noncovalent association of the molecules in the gas phase at 298 K. The calculations of G_solv, G_isom, and G_gas are described below.

Structures

In computing the absolute free energy according to Scheme 1, x-ray structures of the free proteins and their respective complex were used (see Table 1). The ionizable groups in the free proteins and respective complex were protonated or deprotonated according to available experimental pK_a values and the pH of crystallization. All aspartic acid and glutamic acid residues as well as COOH-terminal groups were deprotonated, whereas arginine, lysine, cysteines, and NH₂-terminal groups were protonated. Histidine residues were protonated only if their side-chain nitrogen atoms were within 3.5 Å of a hydrogen acceptor in the crystal structures. However, for the proteins studied here, the histidine side chains were found not to be involved in hydrogen bonding. Thus, they were treated as neutral by protonating at N² or N¹ according to their local environment in the crystal structure. In trypsin, His⁴⁰ and His⁹¹ were protonated at N², but His⁵⁷ was protonated at N¹ (Jackson and Sternberg, 1995). BPTI has no histidine residues. For the Fv region of mouse monoclonal D1.3 antibody and lysozyme, all the histidines were protonated at N² (Tanokura, 1983). Only one histidine is located in the interface region in the complex crystal structures; viz., His⁵⁷ in the trypsin/BPTI complex.

First, hydrogen atoms were added to the crystal structure using the HBUILD module of the CHARMM program (Brooks et al., 1983). The resulting structure was subjected to several steps of minimization using steepest descent followed by adopted-basis Newton-Raphson with strong (10 kcal/mol/Ų) harmonic constraints on all heavy atoms. This relieved close contacts in the protein without disrupting its overall conformation (Philippopoulos and Lim, 1995), as evidenced by RMSDs from the respective crystal structure of 0.064 Å (see Table 2). The relaxed all-hydrogen x-ray structures, denoted by [r]_sln, [l_sln, and [R · L]_sln, were used to compute the solvation free energies (see below). For the isomerization step, each all-hydrogen x-ray structure was energy minimized (initially with constraints, then without) using a dielectric constant of one for 2500 steps of steepest descent, followed by ~7500 steps of adopted-basis Newton-Raphson until the average force was <2 × 10⁶ kcal/mol/Å. The fully minimized [r^min]_gas, [l^min]_gas, and [(R · L)^min]_gas structures, which remain close to their respective crystal structure (see Table 2), were used to obtain the gas-phase complexation energies and entropies. All energy minimizations were carried out using the CHARMM program (Brooks et al., 1983) and the version 22 all-hydrogen forcefield (MacKerell et al., 1998).

In computing G⁰ using the rigid-binding approximation, the starting conformations of the unbound r and l proteins were obtained from the relaxed all-hydrogen x-ray structure of the R · L complex; i.e., [r]_sln = [R]_sln and [l]_sln = [L]_sln. These rigid free structures exhibit heavy atom RMSDs from the respective crystal structures that range from 0.97 to 1.29 Å (see Table 2). The latter value is greater than the motion of protein crystal atoms, 0.3-0.5 Å, derived from their B-factors (Froloff et al., 1997; Luzzati, 1952). Thus, in the rigid-binding approximation, G⁰ is computed according to
The G⁰ was also evaluated using free receptor and free ligand structures that have been modeled from the respective complex by minimizing the receptor R and ligand L conformations using CHARMM with a distance-dependent dielectric constant and constraints on all heavy atoms (denoted by min + cons in Scheme 3). In this case, [r]_sln = [R^min+cons]_sln and [l]_sln = [L^min+cons]_sln. Their heavy-atom RMSDs from the respective crystal structures are similar to those of the rigid structures (see Table 2).

Calculations

Solvation free energy G_solv

(4)

(5)

(6)

(7)

(8)

(9a)

(9b)

Isomerization free energy G_isom

(10)

(11)

Gas-phase binding free energy, G_gas

(12)

(13)

(14)

(15)

(16)

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

In computing G⁰ according to Schemes 1-3, two assumptions were made: the free energy components in Eq. 1 were assumed to be additive; and they were based on a single time-averaged x-ray structure rather than ensemble averaged. Hence, the present scheme is best suited for the noncovalent binding of two proteins that undergo only slight conformational changes to form a high affinity complex, as in the case here (see Table 1). However, for the binding of flexible peptides to proteins, free energy estimates based on a single conformation of the unbound and bound species may not be appropriate.

Tables 3 and 4 summarize the results for D1.3·Hel and trypsin·BPTI binding, respectively. In each case, the binding free energy was computed using three different sets of structures for the free proteins (see Methods): the relaxed all-hydrogen x-ray structures (Scheme 1); rigid [r]_sln = [R]_sln and [l]_sln = [L]_sln structures (Scheme 2); and [r]_sln = [R^min+cons]_sln and [l]_sln = [L^min+cons]_sln (Scheme 3). In what follows, we will first compare the free energies computed using the relaxed all-hydrogen x-ray structures with experiment. Agreement with experiment then allows us to assess the relative contributions from protein-protein versus protein-solvent versus solvent-solvent interactions and electrostatic versus vdW forces to the binding free energy. We then evaluate the accuracy of binding free energies computed using rigid as well as [R^min+cons]_sln and [L^min+cons]_sln free protein structures.

Comparison between computed and experimental G_exp: GB versus PB

Tables 3 and 4 show that the G (11.4 kcal/mol for D1.3·Hel and 18.6 kcal/mol for trypsin·BPTI) based on [r]_sln, [l]_sln, and [R · L]_sln relaxed x-ray structures and G are in excellent agreement with the experimental values (11.4 and 18.1 kcal/mol). In contrast, the G based on [r]_sln, [l]_sln, and [R · L]_sln relaxed x-ray structures and G exhibit larger deviations from the experimental numbers than G, whereas the G and G differ by roughly 20%. The remarkable agreement between G and experiment indicates that systematic errors involved in computing free energies/energies in the right and left legs of Scheme 1 have largely cancelled. Because the G based on x-ray structures are more negative than the respective G by 7.6%, but the CPU time needed to compute G is 5 to 6 times less than that for G, the GB formulation together with Scheme 1 appears to be an efficient and reliable way of computing binding free energies and obtaining trends (see below).

Solvation versus gas-phase contributions to binding affinity

Scheme 1 enables the net binding free energy to be dissected into solvation (protein-solvent and solvent-solvent) versus gas-phase (protein-protein) contributions. The latter, which is a sum of G_gas and G_isom terms, is negative, thus favoring complexation, whereas G_solv is positive and opposes binding (Tables 3 and 4). For both systems, gas-phase vdW and electrostatic protein-protein interactions favor binding (E and E negative). Note that E is negative even though the receptor and ligand have net positive charges (Table 1). The solvent-solvent cavitation term (Eq. 5) also favors protein-protein complexation (G negative) as less work is required to create the complex cavity than the free protein cavities in solution. In contrast to the favorable protein-protein and solvent-solvent interactions, vdW and electrostatic protein-solvent interactions generally oppose binding (G and G positive) due to the cost of desolvating the unbound proteins.

Decomposition of G⁰ into component energies

Tables 3 and 4 show that, for both systems, G (10.4 kcal/mol for D1.3·Hel and 12.0 kcal/mol for trypsin·BPTI) roughly cancels the gas-phase conformational entropy term, TS^conf (9.3 kcal/mol for D1.3·Hel and 11.9 kcal/mol for trypsin·BPTI). This finding can be rationalized if G  TS_solvent, based on the fact that G is related to the hydrophobic effect, which, at room temperature, is dominated by the solvent entropy term. The observed cancellation can then be attributed to the increase in the solvent entropy upon complexation (as protein side chains at the interface release water molecules), and the concomitant decrease in the solute entropy (as protein side chains at the interface lose torsional degrees of freedom upon interacting with other residues) (Novotny et al., 1989). To verify whether the observed cancellation of G and TS^conf is general, these two quantities were also computed for other systems for which x-ray structures are available for both the free and bound proteins (although the structures for these systems may not be as accurate because they do not satisfy all the criteria specified in the Introduction, G and TS^conf are not very sensitive to the structures because they depend on the SASA, see Tables 3 and 4). The results in Table 5 show that G and TS^conf oppose each other and their difference is less than 1.5 kcal/mol, indicating that G and TS^conf generally offset one another. It should be emphasized that the observed G and TS^conf compensation rests on the magnitude of the coefficient , 7.2 cal/mol/Ų, used in this work (see below).

View this table: [in this window] [in a new window]   TABLE 5   Changes in the side-chain conformational entropy (TSconf) and nonelectrostatic solvation free energy (G) in various protein-protein associations

Because the G and TS^conf terms offset one another in Tables 3 and 4, the remaining components of the binding affinity can be partitioned into net gas-phase protein-protein vdW effects (E^vdw), net protein-protein and protein-solvent electrostatic effects (G^elec), and the sum of translational, rotational, and vibrational entropy changes (TS^{trans+rot+vib}). The magnitudes of E^vdW, G^elec, and TS^{trans+rot+vib} in Tables 3 and 4 are sizable, indicating that all three components contribute to the overall binding affinity, so that

(17)

The G computed using Eq. 17 are 9.7 kcal/mol for D1.3·Hel and 17.9 kcal/mol for trypsin·BPTI, which underestimate the experimental numbers by 1.7 and 0.2 kcal/mol, respectively. Of the three terms in Eq. 17, only the gas-phase protein-protein vdW interactions favor binding (E^vdw negative). Furthermore, the magnitude of E^vdw is larger than that of G^elec or TS^{trans+rot+vib} for both D1.3·Hel and trypsin·BPTI complexation. These findings are consistent with the experimental observation that enthalpy drives D1.3·Hel binding from measurements of enthalpy and entropy changes by titration calorimetry (Bhat et al., 1994). Hence, close packing and shape complementarity play important roles in noncovalent protein association. The net electrostatic interactions generally oppose formation of complexes with an interface area between 1400 and 1700 Å² (G^elec positive). This finding does not mean that electrostatic interactions are unimportant for protein binding, but shows that the gain in hydrogen bonding and charge-charge protein-protein interactions across the interface is offset by the loss of favorable electrostatic protein-solvent interactions. Although translational and rotational entropy loss inevitably oppose complexation (S^trans+rot negative), vibrational entropy can favor or disfavor protein-protein association: S_vib is negative for D1.3·Hel complexation, but is positive for trypsin-BPTI binding (see also Discussion).

Because the magnitude of the TS^trans+rot term is similar for two ligands of roughly equal volumes, l and l', binding to a common receptor r, their binding free energy difference will be given by

(18)

where the prime denotes the thermodynamic change for binding of l'. Hence, the discrimination between two similarly shaped ligands, which defines specificity, is governed by close packing of the protein surfaces so that the proper hydrogen bonds can be formed (McCammon, 1998). Eq 18 shows that the affinity of a drug ligand l' for its receptor protein with an interface area between 1400 and 1700 Å² can be improved by mutations that increase the magnitude of E'^vdw, but decrease the magnitude of G'^elec.

Evaluation of the rigid-binding approximation

The RMSDs of the heavy atoms in the free x-ray structure from the corresponding complex crystal structure is <1.3 Å for both systems studied (Table 1). Therefore, the structural changes upon complexation are relatively small and it seems reasonable, to a first approximation, to neglect any conformational changes upon protein-protein association; i.e., to assume [r]_sln = [R]_sln and [l]_sln = [L]_sln. In fact, previous studies used this rigid-binding approximation to compute the free energy for D1.3·Hel and trypsin·BPTI binding (Tables 3 and 4). The binding free energy computed using Scheme 2 does not agree with the experimental value (Tables 3 and 4, column 3), and is even quite positive for the binding of D1.3 antibody to Hel. The discrepancy between theory and experiment arises from both protein-protein as well as protein-solvent electrostatic interactions. The rigid free receptor and free ligand conformations isomerize to local minima, which have less favorable electrostatic interactions than the local minima derived from the relaxed x-ray structures of the free proteins. Consequently, the E values computed using the rigid free receptor and free ligand structures (242 and 315 kcal/mol) are more negative than the respective numbers in the second column of Tables 3 and 4 (106 and 228 kcal/mol). In contrast, the G values computed using the rigid free structures (119 and 198 kcal/mol) are more positive than the respective numbers in column 2 of Tables 3 and 4 (98 and 85 kcal/mol).

In the case of D1.3 antibody binding to Hel, the relative component contributions to the binding free energy computed using Scheme 2 do not agree with those obtained using Scheme 1. In particular, Scheme 2 predicts that the dominant contribution to the net binding free energy for D1.3·Hel complexation is not E^vdW, but the entropy term, TS^{trans+rot+vib}, in contrast to Scheme 1 and experimental observations (see above). These results show that, even if the structural changes upon complexation are known to be small, the rigid-binding approximation (Scheme 2) will generally not yield accurate binding free energy, and may also not reveal the dominant contribution to the binding free energy.

Evaluation of free R^min+cons receptor and free L^min+cons ligand structures

In cases where the free receptor or free ligand structures have not been experimentally solved but their complex crystal structure is known, the former can be predicted from the latter by minimizing the receptor and ligand conformations in the complex (see Structures section above) provided that binding results in minimal conformational changes. Like the rigid structures, the [r]_sln = [R^min+cons]_sln and [l]_sln = [L^min+cons]_sln structures predict a positive binding free energy for D1.3·Hel complexation, whereas they severely overestimate the magnitude of the free energy for trypsin·BPTI binding. However, they yield trends in the relative component contributions to the binding free energy that are similar to those based on the relaxed x-ray structures (Tables 3 and 4). The magnitude of the component contributions based on the [R^min+cons]_sln and [L^min+cons]_sln structures are in-between the respective numbers based on the x-ray and rigid structures (Tables 3 and 4). This suggests that approximating [r]_sln and [l]_sln by [R^min+cons]_sln and [L^min+cons]_sln, respectively, in cases where slight conformational rearrangements accompany complexation, as in the cases studied here, can reveal qualitative features in the binding free energies.

Dependence of component free energies on free and complex structures

The results using Schemes 1-3 show that accurate structures are needed to obtain reliable absolute binding free energies (Sharp, 1998). This sensitivity of the atomic coordinates is expected because electrostatic and vdW interactions depend on the interatomic distances. Hence, the E^vdw and G^elec terms computed using different free structures vary significantly (see Tables 3 and 4). The vibrational entropy term is also sensitive to the local minimum structure used to compute the frequencies (see Table 6). In contrast, the G and TS^conf terms computed using different free structures do not vary as much as they depend on the solvent-accessible solvent surface area. These findings are in accord with previous works (Froloff et al., 1997; Jackson and Sternberg, 1995; Novotny et al., 1997).

View this table: [in this window] [in a new window]   TABLE 6   Lowest four nonzero vibrational frequencies (cm1) and vibrational entropies (TS in kcal/mol) from normal mode calculations for the fully minimized structures*

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Below, we compare our above findings with results obtained in previous studies for the same or homologous systems and for other protein-protein complexes.

Comparison of the various free energy decomposition schemes

The scheme used here to compute the free energy for protein-protein complexation is closest in spirit to that used by Jayaram et al. (1999) for protein-DNA complexation. The key difference lies in the isomerization step in Scheme 1, which was introduced to obtain local-minimum structures (with no imaginary frequencies) for normal mode frequency calculations. In contrast, Jayaram et al. computed first the isomerization/adaptation of the protein structure in the free state to the conformation in the complex state (the adaptation energy), and, subsequently, their desolvation (see Scheme 4). As a result, the gas-phase complexation step corresponds to species that are not in energy minima, thus their normal mode frequencies cannot be computed.