INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

The interface between two proteins in a complex is generally as tightly packed as the protein interior (LoConte et al., 1999). Direct interactions and the removal of water from the interface provides the attractive contributions to the binding free energy that more than compensate for the loss of translational, rotational, and side chain entropy upon association. In addition, for protein pairs such as barnase and barstar, the highly specific electrostatic interactions provide long-range steering effects, resulting in an association rate that is higher than it would be without such attractive forces (Schreiber and Fersht, 1996). We have recently shown that desolvation also provides specific attractive interactions that, albeit shorter-range and weaker than the electrostatic steering, can also increase association rates by several orders of magnitude (Camacho et al., 1999, 2000b).

The tightly packed interface and its thermodynamic and kinetic consequence are all lost when a complex is constructed from the separately determined structures of two proteins. Indeed, more often than not one finds some interfacial structural motifs in "wrong" positions (Koshland, 1958, 1963, 1994; Jorgensen, 1991), resulting in steric clashes and unfavorable electrostatic interactions, even for high resolution x-ray structures and for proteins whose backbone remains practically invariant in the process of binding (Jackson et al., 1998; Vakser et al., 1999; Camacho et al., 2000a). These findings emphasize that upon binding, the protein interactions should lead to some degree of induced fit (Koshland, 1994; Jorgensen, 1991), resulting in the tightly packed interface. Induced fit theory as introduced by Koshland (1958) explains the origin of functional specificity and describes how a substrate changes the structure of an enzyme to bring its catalytic groups into the proper alignment, whereas a nonsubstrate does not. However, the notion of induced fit is also frequently used in a more general sense as the origin of binding specificity, i.e., the collection of conformational changes resulting in optimal interactions when two molecules come in contact (Jorgensen, 1991).

Though induced fit is necessary for forming the well-packed interface seen in protein-protein association, it cannot fully explain important kinetic requirements of the binding process. For instance, if side chains are in wrong conformations before the two molecules contact each other, then the long and medium-range forces due to electrostatics and desolvation that normally bring the two molecules together can be completely eliminated (see below). Furthermore, if these misfolded side chains should rearrange in order to restore the proteins attraction, then these changes must occur within the time scale of an encounter between the receptor and ligand. The lack of nonspecific protein aggregation constrains this time scale to a few nanoseconds (Northrup and Erickson, 1992; Camacho et al., 2000b). On the other hand, even a small rotation of a medium-sized side chain in the protein interior can take as much as 1 s (Creighton, 1993). In this context neither the traditional lock and key (Fischer, 1894) nor the induced fit (Jorgensen, 1991) model seems to provide a rationalization for the attraction between proteins required for a fast and specific binding process. In particular, it becomes an open problem what positions the key side chains occupy before and during the process of protein-protein association.

We use a novel molecular dynamics (MD) algorithm to analyze the coupling between water molecules and flexible surface side chains in three protein ligands. MD has been used to analyze the short time scale behavior of small molecules, mainly water, around typically rigid macromolecules (Brunne et al., 1993; Steinbach and Brooks, 1993; Makarov et al., 1998). However, it has also been observed experimentally that water-protein interactions influence protein dynamics (Zanotti et al., 1999). Recognizing that for many protein-protein complexes, surface side chains are the basis for the detailed chemistry that gives rise to their fast and highly specific association, we studied the effects of explicit solvent on the conformational dynamics of surface side chains from the independently solved (unbound) structures of barstar, bovine pancreatic trypsin inhibitor (BPTI), and lysozyme in the absence of the receptor using MD simulations.

This paper revisits the principle, proposed nearly a hundred years ago by Langley (1907), that a ligand is attracted to its receptor by a specific pre-existing affinity between the molecules. In particular, we considered barstar, BPTI, and lysozyme because they show little or no affinity for their receptors when calculated using the x-ray or nuclear magnetic resonance (NMR) structures of the unbound ligands. We focused on the side chains found at the binding interface of the complex structures. Our simulations show that by properly solvating these side chains, the rotamers sampled are frequently those seen in the complex. We also checked that the use of these rotamers restores the expected affinity between the molecules. This confirms that for such side chains, a localized lock-and-key-like complementarity exists before any interdigitation of the monomers. In contrast, side chains close to the perimeter of the interface, though important for the stability of the complex, are generally found in nonspecific rotamers, suggesting that their final conformations are induced by the presence of the receptor.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

Protein complexes

We chose to study representative systems from each of the three major classes of protein-protein complexes for which structures are known. These are the barnase-barstar complex (PDB code 1brs) from the RNase-inhibitor family, the trypsin-BPTI (2ptc) and the trypsin-kallikrein (2kai) complexes from the protease-inhibitor family, and hen egg-white lysozyme bound to two different antibodies, Fab D44.1 (IgG1, ) (1mlc) and HyHEL5 (2hfl) from the antigen-antibody family (Fig. 1). The structures of these receptors and ligands have also been solved independently and, hence, have been analyzed in numerous docking studies (Jackson et al., 1998; Camacho et al., 2000a; Vakser and Aflalo, 1994; Vakser et al., 1999).

View larger version (42K): [in this window] [in a new window]   FIGURE 1   Receptor/ligand interfaces in the complex structure with the simulated side chains. In blue are the side chains found in the monomer (unbound) ligand structure. In red are side chains found in the cocrystallized (bound) complex. In yellow are nearby charged residues on the cocrystallized receptor. Left: the barnase-barstar complex. Note that Glu76 of 1bta corresponds to Glu74 of 1brs. Middle: the trypsin-BPTI complex. Right: the Fab arm of antibody IgG1 D44.1, , complexed with hen egg-white lysozyme.

MD simulations were carried out on the unbound protein ligands alone. We focused our attention on two charged side chains on each ligand, selected on the basis of their involvement in interactions with the receptor. The first of the side chains in each structure is fully buried in the complex ("interfacial" in Table 1), whereas the second side chain is always partially exposed to the solvent ("peripheral"). In order to improve the efficiency of the simulations, all other side chains were kept fixed in their unbound conformations (see Table 1 and Fig. 1). We have checked that our results and conclusions did not significantly change when we allowed neighboring side chains to move (see MD protocol). As shown below, the conformations of the selected side chains were substantially affected by the water-ligand interactions. In contrast, the interfacial side chains on the unbound receptors in our examples (e.g., 1bao, 2ptn, 1mlb) differed by <1 Å root mean square deviation (RMSD) between their free and bound states.

Surface of active polarons: semi-explicit solvation model

The MD simulations presented in this study were carried out using a semi-explicit solvent model formulated and developed in previous papers (Brower and Kimura, 1998; Kimura et al., 2000). Briefly, the method involves using a layer of explicit water surrounding the solute and a set of surface charges that are placed directly on the oxygen atoms of the explicit water molecules. The cooperative effects of these charges recreate the polarization of the external bulk not included in the simulation. Here, the external water is modeled as a dielectric continuum with ₁ = 78.0. The surface charges are updated continually throughout the simulation using a discretized version of the self-consistent relation,

(1)

This expression is the same as that used in the boundary element method (Zauhar and Morgan, 1985, 1990; Purisima and Nilar, 1995) and is derived from Poisson's equation and appropriate boundary conditions (Kimura et al., 2000). In the equation, ₀ is the dielectric constant of the vacuum, r_k and q_k are the coordinates and magnitudes of the explicit charges of the system (including both solute and explicit water),  are the positions of the oxygens of the explicit water molecules, are the outward normal vectors on these oxygens, and the integration is performed numerically around a closed boundary that is approximated by the oxygen positions of the layer of water surrounding the solute of interest. The first term in this equation accounts for the potential due to the explicit charges in the system, and the second comes from the potential due to the surface charge distribution on the boundary. We refer to these water molecules with their extra charge as a Surface of Active Polarons (SOAP). In addition to the extra charges, short-range forces and thermal fluctuations are applied to the SOAP particles to model hydrostatic pressure and transmission of energy from the bulk. The model has been shown to reproduce solvation energies of 9 ions and 14 polar and charged amino acids with high accuracy (correlation coefficients with experiment 1.000 and 0.995, respectively; Kimura et al., 2000).

Initial conditions

The structures for barstar, BPTI, and lysozyme were obtained from the Protein Data Bank (PDB codes 1bta, 6pti, and 1lza) and modeled using the OPLS/AMBER united atom force field (Jorgensen and Tirado-Rives, 1988) supplied with the TINKER (Pappu et al., 1998) v3.6 distribution. The molecules were solvated with a 9.0-Å layer of explicit TIP3P (Jorgensen et al., 1983) water and minimized until the gradient of potential reached 5.0 Kcal/mol/Å with all solute atoms except the side chains of interest fixed. A very brief initial MD simulation was carried out (4000 steps for 1bta, 3000 steps for 6pti, and 5000 steps for 1lza using a 2-fs time step) to obtain the initial surface charges on the explicit water. Data from these preliminary simulations were not used in our analysis for this paper.

MD protocol

All MD simulations were carried out using the TINKER package (Pappu et al., 1998), which was modified to include the semi-explicit solvation model described above and also a newly parallelized force loop. The modified package was compiled and executed on an Origin 2000 (SGI, Mountain View, CA). The complete coordinate set was saved every 100 fs during each simulation for off-line processing. Water bond lengths were constrained using the rattle algorithm (Andersen, 1983), and all solute atoms other than the two side chains of interest were rigidly constrained to the original coordinates of the unbound ligand PDB structures. Coordinates were updated using the velocity Verlet algorithm with a 2-fs time step. For the electrostatics, a 20 Å distance cutoff was used with cubic spline potential smoothing to avoid potential discontinuities.

To check that fixing neighboring residues did not adversely affect our results, we repeated the lysozyme simulation for 500 ps, allowing the residues Arg45, Asn46, Thr47, Asp48 and Gly67, Arg68, Thr69, Pro70 to move while keeping everything else fixed (not shown). Conformational sampling of the two side chains of interest for this system, Arg45 and Arg68, showed slower convergence but very similar trends, as when the neighboring residues were fixed.

Grid potential

We optimized the simulations by precalculating the effects of fixed atoms under a plane located 6.0 Å below the lowest atom of the simulated side chains. The electrostatic potential was computed at points sampled on a cubic lattice with 0.5 Å spacing within an approximately 60 × 50 × 60 Å box that enclosed the part of the molecule above the plane. This potential was computed using both the solute atoms and the solvent particles with their extra charge, which was calculated and saved during the short initial simulation. This ensured that the effects of the dielectric boundary at the far end of the molecule was approximated in the grid potential as a static reaction field due to the charges of the system. The resulting grid potential was used in all MD simulations using linear interpolation to approximate the continuously varying potential between grid points. In addition, a reflective barrier was placed at the plane separating the implicit and explicit parts of the molecule to prevent water from escaping below the plane. This simplification improves the efficiency of the simulations by approximately a factor of 9.

Evaluation of binding free energy

As in a recent analysis by Lee et al. (2000), we used an effective free energy function to estimate energies for conformations extracted from our MD simulation to compare with that of the crystal complex. Specifically, we computed the binding energies for the cocrystallized protein structures, the complex formed by the native receptor and unbound (independently solved) ligand, as well as for solvated side chain conformations identified in our MD simulations.

The binding free energy is calculated by the expression (Vajda et al., 1994; Novotny et al., 1989)

(2)

where E_coul denotes the direct electrostatic energy between receptor and ligand, G_des corresponds to the desolvation energy including side chain entropy loss, and G_rottrans accounts for translational, rotational, vibrational, and cratic effects. E_coul is calculated using a distance-dependent dielectric (McCammon et al., 1979) equal to 4r, enforcing a minimum atom-to-atom distance separation equal to the sum of their corresponding van der Waals radii to avoid artificial overlaps. G_des is estimated using an atomic contact energy (ACE; Zhang et al., 1997). For the most part, G_rottrans accounts for the loss of rotational and translational degrees of freedom upon binding. This entropic barrier opposing protein association must be surmounted in order for proteins to bind. G_rottrans is typically assumed to be a weak function of the size and shape of the interacting proteins (Vajda et al., 1994; Novotny et al., 1989). Here, we assume G_rottrans = 5.0 Kcal/mol, a value that fits our data well.

Eq. 2 will be used to understand how the side chain conformations affect the binding free energy. In all calculations the receptor is kept in its bound structure, but the conformation of the ligand varies. For example, G_bind^calc(ul) and G_bind^calc(bl) will denote the calculated free energies of the bound (b) and the unbound (u) conformations of the ligand (l). Furthermore, within the rigid-body approximation, the apparent loss of affinity G_loss for the unbound ligand is estimated as

(3)

The same equation (Eq. 2) for estimating changes in binding free energy will be used when key side chains on the ligands are replaced by conformations identified through MD simulation (see Results and Table 3).

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

Free energy calculations using x-ray structures

Table 2 shows a list of the calculated and observed binding free energies, G_bind^calc and G_bind^exp, respectively. We also list G_loss defined in the methods, and G_loss^fix, which denotes the loss of affinity of the unbound ligand when the two selected side chains are replaced by their conformations found in the crystal structure of the complex. The agreement between G_bind^calc and G_bind^exp is good, with the single exception of the lysozyme-HyHel5 complex (3hfl/1lza), which other groups have also found hard to model (Novotny et al., 1989). Including the water molecules found at the interface of the antibody-antigen complexes yields an improved estimate of G_bind^calc = 17 Kcal/mol for 3hfl/1lza (unpublished). This, together with previous validation of ACE (Zhang et al., 1997), gives credence to our estimates of G_loss.

Because the side chains selected for our study are charged and form salt bridges in the complex, the binding affinity is particularly sensitive to their position. We find that the loss of affinity amounts to >75% of the total binding free energy of the complex. This loss is significantly reduced when the two selected side chains are in the rotamer conformation found in the complex, while all the other side chains remain in their unbound conformations (see G_loss^fix in Table 2). It is important to emphasize that error bars on G_loss do not change the main observation that misfolded side chains can significantly affect the receptor/ligand binding affinity, both for the complex and for near-native structures (Camacho et al., 2000a).

Dynamics

From the analysis of the nanosecond MD simulations of the solvated ligands shown in Fig. 2, we conclude that the side chains fully buried in the complex structure (Asp39 for 1bta, Lys15 for 6pti, Arg68 for 1lza) are frequently in a state within 1.5 Å RMSD of the side chain conformation in the complex structure. In contrast, the side chains that are partially exposed in the complex (Glu76 for 1bta, Arg17 for 6pti, Arg45 for 1lza) do not sample their native-like complex conformation. To quantify the sampling of the different rotamers, each side chain conformation was clustered around Dunbrack's rotamer library (Dunbrack and Cohen, 1997) and ranked in Table 3.

View larger version (34K): [in this window] [in a new window]   FIGURE 2   Heavy atom RMSD of each side chain studied relative to the corresponding side chains found in the complex. RMSD values are plotted every 50 steps during the MD simulation. The dashed line shows the RMSD of the initial unbound monomer structure.

Buried side chains move to positions favorable for protein-protein recognition

Strikingly, Asp39 of 1bta (Fig. 2, top) very rapidly moves and settles into a conformation close to that of the complex (approximately 1.2 Å RMSD), in the absence of the interacting receptor. The simulations show only one dominant cluster (Table 3) containing the native-like rotamer (no. 4) of Asp 39 in the complex crystal structure (1brs). This side chain cluster alone, in the unbound ligand PDB structure, yields a substantial improvement in ligand/bound-receptor interaction energy, described by the mean gain G_gain^Asp39 = G_loss^Asp39G_loss = 7.1 Kcal/mol (see Table 2), when compared to the binding affinity of the intact unbound ligand and the bound receptor. In this expression G denotes averaging the free energy over the cluster.

Although we see more fluctuations for Lys15 of BPTI, this side chain also frequently visits essentially the conformation in the complex at 0.7 Å RMSD. Out of all the possible conformations, the rotamers in the unbound and bound structures rank first and second in their frequency of sampling. Consistent with the interaction energy in the complex, the cluster of structures closer to the native side chain have an energy improvement of G_gain^Lys15 = 8.1 Kcal/mol.

For lysozyme, the unbound crystal structure shows Arg68 and Arg45 in electrostatically favorable pockets on the protein surface. Despite its large fluctuations between 1 and 6 Å RMSD, upon solvation Arg68 repeatedly returns to a state very close to the complex (rotamer no. 14) near 1.0 Å RMSD and spends as much as 20% of the time in conformations with low (<2.0 Å) RMSD values. Examination of these clusters reveals that they are structurally very close to the bound conformation, and a small adjustment of about 1 Å (perhaps caused by the approaching antibody) is all that is required to bring this side chain to its final bound position. Although this side chain has an initial unbound monomer conformation that is already close to the complex (Fig. 1), the fact that it returns to these low-RMSD conformations suggests the existence of an attractive well in that region that competes favorably with the solvation entropy. Because the unbound rotamer and the rotamer found in the native complex state are very similar, there is no significant energy improvement for this side chain in its solvated conformation.

Partially solvent exposed side chains are typically found in nonspecific rotamer conformations

For Glu76 of 1bta, there is a dominant cluster somewhat close to the rotamer found in the unbound structure (rotamer 27), and only a very brief visit to a rotamer close to the one found in the bound structure (rotamer 22).

The largest cluster for Arg17 of BPTI is near the rotamer in the initial unbound structure; we do not see sampling of the native complex rotamer. Nevertheless, we find substantial energy improvements for many of the nonspecific highly ranked clusters, suggesting that the rotamer in the unbound crystal is in a particularly unfavorable conformation.

Finally, Arg45 in 1lza appears to make only rare visits to a 3.5 Å RMSD baseline, and otherwise stays further away from the initial monomer structure, near 5.0 Å RMSD. Arg45 neither includes the bound (rotamer 33) or unbound (rotamer 59) rotamers in its frequently visited clusters and yet we see improvement in energies for the top two clusters of about G_gain^Arg45 = 7 Kcal/mol. These clusters are dominated by solvation effects, suggesting that the rotamer found in the unbound structure is unfavorable in solution. The main observation here is that in the original unbound structure, Arg45 completely blocks the binding site, preventing any possible affinity between the monomers. Our MD simulations suggests that in solution, this side chain moves away from the binding area, allowing the proper complementary surfaces to interact.

	DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION REFERENCES

The goal of this study was to investigate the conformational dynamics of side chains on protein surfaces that play an important role in complex formation (for an animated visualization see http://engpub1.bu.edu/~srk/keys.html). The hypothesis is that key side chains must already be in positions suitable for binding in the absence of their receptors to allow for the recognition process to take place in a biologically feasible time scale, and to yield enough intermolecular affinity to overcome the entropic barrier to binding. Among the side chains we studied, the three that are located deep within the binding interface indeed exhibit key-like behavior. For Asp39 of barstar, Lys15 of BPTI, and Arg68 of lysozyme, we find highly populated clusters of conformations that are within 1.5 to 2 Å RMSD from the structures found in the complex. This sampling occurs within the nanosecond time scale, which also corresponds to the typical lifetime of an encounter complex (Northrup and Erickson, 1992; Camacho et al., 2000b). Thus, we conclude that in all likelihood when ligands interact with their receptors their key side chains are in their "right" bound-like conformations. Table 2 shows that the interactions of the appropriately oriented side chains with their substrate significantly enhances the binding affinity between monomers. The highly specific nature of the interactions between these key residues and their receptors also suggest that mutations of these side chains should likely result in a decrease of the association rate. The latter is consistent with the observation (Castro and Anderson, 1996) that the Lys-15-Ala mutation on BPTI leads to a 200-fold decrease of the on rate of the BPTI-trypsin association. From a thermodynamic point of view, we have also checked that an alanine mutation on any of the side chains considered here results in an overall decrease of the binding affinity of the complex.

For the other three side chains on the periphery of the binding interface, we find highly populated nonspecific conformations that differ from the initial unbound structures. In the original unbound ligands, these side chains are often found interfering or even blocking the binding site; we find that solvation effects remove these chains from the binding area to positions accessible to the final complexed state. Our simulations suggest that these side chains act more as "latches" that hold the molecules together rather than as keys fitting snugly into the receptor. Consistent with an induced fit mechanism, these latches should be open before docking and only fasten on to their salt-bridge partners late in the binding process. An implication of these results is that these salt bridges should mostly yield a decrease in the off-rate as opposed to an increase in the association rate. The latter is consistent with the mutagenesis experiments of Schreiber and Fersht (1996) where the mutation of Glu76 to alanine in barstar resulted in a very modest decrease of the association rate with barnase. Similarly, Castro and Anderson (1996) have shown that the mutation of the latch residue Arg17 in BPTI to alanine again leaves the on rate almost unchanged, whereas the off rate is larger than for the wild-type.

Crystallographic data shown in Fig. 3 provide further support for the generality of our findings. Namely, the key-like feature of the buried side chains is confirmed by the fact that BPTI complexed with kallikrein has the same Lys15 rotamer as BPTI with trypsin. The same is found for Arg68 in lysozyme complexed with antibodies D44.1 and HyHEL5. The latter is true despite the fact that these complexes have very different binding free energies and association rate constants. On the other hand, the induced fit mechanism suggested for the nonspecific partially exposed side chains would suggest that their final complex conformations should not necessarily be the same for different complexes. This is indeed the case for the bound rotamer of Arg17 in BPTI complexed with trypsin and kallikrein, and for Arg45 in lysozyme complexed with antibodies D44.1 and HyHEL5 which, as shown in Fig. 3, are found in very different rotamers.

View larger version (30K): [in this window] [in a new window]   FIGURE 3   Interfacial side chain conformations of BPTI and lysozyme when bound to two different receptors. For BPTI (top row), in blue and red are side chain conformations found in the bound complex with trypsin (2ptc) and kallikrein (2kai), respectively. For lysozyme (bottom row), in blue and red are the conformations found in the complex with antibodies D44.1 (1 mlc) and HyHel5 (3hfl), respectively.

The binding mechanism that emerges from our analysis is one in which the chemical affinity between receptors and their ligands is crucial for a fast recognition process. This affinity is enhanced when key side chains, most important for binding, are properly oriented near their conformations seen in the complex, allowing for the rapid and highly specific interdigitation of the molecules. This mechanism is consistent with the thermodynamic maps (Camacho et al., 1999) of electrostatic and desolvation free energies of receptor/ligand systems, which clearly show a funnel-like shape around a properly aligned precursor state. Furthermore, Brownian dynamic simulations (Camacho et al., 2000b) have confirmed that these molecules could efficiently avoid a lengthy interdigitation process over large contact areas by locking into a well-defined diffusion accessible precursor structure. Partially solvent exposed side chains do not appear to play this key-like role but they contribute to the stability of the complex.

An interesting practical application of our finding is to use the identification of conformational preferences of surface side chains in combination with standard docking algorithms to predict the structure of the complex. Challenges that remain for this type of procedure include the reduction of the number of side chain conformation candidates to a reasonably small number and the inherent computational cost of atomistic MD or Monte Carlo simulation using an accurate solvent model. The benefits gained from this approach are evident, as it decouples the docking procedure into separate side chain and rigid-body searches, thereby reducing the space of possibilities to a manageable size.

In summary, we find that important information can be learned from the characteristic distribution of rotamers resulting from the coupled dynamics of flexible surface side chains and its hydration layer. Our findings suggest that the unique environment of the side chains consisting of neighboring atoms and the solvent, which is ultimately encoded in the primary structure, holds the keys for protein binding. Complex-forming proteins appear to have interfacial side chains that frequently sample favorable conformations in order to meet the stringent kinetic requirements for binding.