| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
1 Department of Physics and
2 Department of Biochemistry, University of Washington, Seattle, Washington 98195, USA
Reprint requests to: David Baker, Department of Biochemistry, University of Washington, Seattle, WA 98195, USA; e-mail: dabaker{at}u.washington.edu; fax: (206) 685-1792.
(RECEIVED November 18, 2003; FINAL REVISION March 2, 2004; ACCEPTED March 2, 2004)
 |
Abstract |
|---|
 TOP Abstract IntroductionResultsDiscussionReferences |
|---|
Keywords: protein–protein interactions; diffusion-limited association rates; orientational constraints; rotational diffusion; long-range interactions; Brownian dynamics
Article published online ahead of print. Article and publication date are at http://www.proteinscience.org/cgi/doi/10.1110/ps.03517304.
|
Introduction |
|---|
|
TOPAbstractIntroductionResultsDiscussionReferences |
|---|
DR (where D is the relative translational diffusion constant and R is the sum of the radii), which yields rates of 109–1010 M–1 sec–1 for associations relevant to proteins. Usually, however, proteins exhibit a highly anisotropic distribution of reactivity over their surface. This can be modeled by localized reactive sites on the surface of the proteins that have to be sufficiently precisely aligned for the complex formation to occur. Purely probabilistic models have tried to account for such steric constraints by multiplying the Smoluchowski rate for uniform spheres by the probability that, in a random encounter, the two molecules are properly aligned ("geometric rate"; Janin 1997 gives an example of this method). This yields rate constants that are typically several orders of magnitude lower than the Smoluchowski diffusion-limited rate and are usually much smaller than the values experimentally observed for biological complexes. It has been found that this discrepancy can be moderated by taking into account the effect of rotational diffusion (Shoup et al. 1981; Northrup and Erickson 1992); additional rate enhancements are brought about by the presence of attractive interparticle forces ("electrostatic steering"; see Schreiber and Fersht 1996; Gabdoulline and Wade 1997; Vijayakumar et al. 1998), and the formation of a weakly specific, loosely bound encounter complex that subsequently evolves into the final bound state (Selzer and Schreiber 1999; Camacho et al. 2000).
To replace the estimation of protein–protein association rates via the geometric rate by a more accurate method, most authors have pursued a computational approach by carrying out explicit numerical simulations of the diffusional association of macromolecules, commonly referred to as Brownian dynamics (BD) simulations (for an excellent review, see Gabdoulline and Wade 2002). Here, the protein molecules are modeled in varying detail, from a simple spherical approximation up to full atomic detail. In the simulation, the molecules are initially placed in random orientations at a fixed initial separation b. Diffusional trajectories, with or without the presence of an interparticle force (such as electrostatic interactions), are then generated by means of the Ermak–McCammon algorithm (Ermak and McCammon 1978). A trajectory is ended either when the molecules have come together in proper orientation to successfully form a complex, or when their separation has exceeded a certain truncation value c > b such that the probability for an encounter has become vanishingly small. The fraction of "successful" trajectories is then used to compute the association rate kon.
Northrup and Erickson (1992) have used such a BD simulation to compute the association rate of spherical molecules with a reactive patch, consisting of four contact points in a 17 Å x 17 Å square arrangement on a plane tangential to the surface of the molecules. A reaction is then assumed to occur if three of the four contact points are correctly matched and within a specified maximum distance. In the absence of any interparticle forces, the authors find an association rate of kon = 105 M–1 sec–1, about two orders of magnitude higher than the geometric rate. Gabdoulline and Wade (2001) compute association rates for five protein–protein complexes using full-atom structures in the presence of long-range electrostatic forces. The reaction condition is defined by formation of subsets of the polar contacts observed in the native complex structure.
In this paper, we present a different route toward estimation of rates of bimolecular association. Instead of using a computer-simulation-based approach such as the method of BD simulations outlined above, we use a recently derived analytical expression (Schlosshauer and Baker 2002) for the association rate of two spherical molecules with anisotropic reactivity in the absence of any interaction forces. The reaction condition is formulated by specifying the ranges of mutual orientations of the two molecules for which complex formation will occur. We thus do not require an exact mutual alignment of the binding partners, but instead assume that favorable short-range interactions "guide" the molecules into their final bound configurations once the molecules are oriented within specified angular tolerances (see Fig. 1
). These tolerances can therefore be viewed as an implicit modeling of attractive short-range forces. We derive estimates for the tolerances from free energy landscapes obtained by sampling configurations within and surrounding the native binding funnel. These values are then used in our analytical expression to compute the corresponding association rates. By determining the size and geometry of the aperture in phase space that must be entered for binding to occur, and rigorously solving the problem of diffusion through this aperture, our approach provides a physically transparent complement to BD simulations for computing binding rates from structures of protein–protein complexes.
|
|
Results |
|---|
|
TOPAbstractIntroductionResultsDiscussionReferences |
|---|
, three ingredients are required. The first is a general theory for computing the diffusion-limited association rate as a function of the orientational constraints associated with properly aligning the two binding sites (the size and shape of the reactive zone in Fig. 1). The second is a method for transforming a diffusion in a potential problem into a free diffusion problem—in the context of Figure 1, an estimate of how deeply the reactive zone lies within the binding funnel (i.e., how far molecules must descend into the binding funnel before they are effectively captured). The third is a method for mapping the binding funnel for two proteins given the crystal structure of the protein–protein complex. We address these issues in the next three sections and then, in the fourth section, use the results to compute approximate diffusion limited association rates from the structures of protein–protein complexes.
Theory for the diffusion-limited association rate with general orientational constraints
Here, we restrict ourselves to a brief outline; the full derivation of our expression for the association rate constant in the presence of general orientational constraints can be found in Schlosshauer and Baker (2002).
We consider translational and rotational diffusional motion of two spherical molecules A nd B with radii RA and RB. To derive an expression for the association rate constant, we solve the steady-state translational–rotational diffusion equation describing the diffusional motion of the two spheres, subject to a reaction condition that ensures that binding can only occur if the mutual orientation of the two spheres is sufficiently close to the orientation in the bound configuration that defines the optimal alignment.
The reaction condition is implemented as follows (see Fig. 2): The centers of "reactive patches" on the two spheres are defined by the intersection of the center-to-center vector with the surfaces of the spheres in the native bound configuration. Each sphere carries its own body-fixed coordinate system {xs, ys, zs}, s = A, B, where the zs axis points at the center of the reactive patch. The angles 
A and B then quantify the distance of the center of each reactive patch to the center-to-center vector, whereas the Euler angles 
and 
denote relative torsional angles between the body-fixed coordinate systems (xA, yA, zA) and (xB, yB, zB).
|
A in addition to the polar angle A to fix the location of the reactive patch on the surface of sphere A. For the formulation of the reaction condition, however, four angles suffice, because the position of the center of the reactive patch is automatically specified through its coincidence with the zA axis. This leaves only one free parameter, namely, the "width" of the patch, which is described by the angle A.
The optimal alignment is then defined by A = B = = = 0 (additionally, the length r of the center-to-center vector must be equal to the sum of the radii of the spheres). Our reaction condition requires that all these angles are sufficiently close to zero for the reaction to occur; that is, the following conditions are fulfilled:
 | (1) |
Using the constant-flux approximation introduced by Shoup et al. (1981), we obtain for the association rate constant (Schlosshauer and Baker 2002):
|
| (2) |
where D = DtransA + DtransB is the (relative) translational diffusion constant, a0 = (4)300(1 – cos 0A)(1 – cos 0B), qllAlB = (2l + 1)(2lA + 1)(2lB + 1)/163, and 
quantifies the extent of diffusion control in the reaction. Furthermore,
 |
where dlmn() denotes the Wigner rotation function. 
is the Wigner 3-j symbol, and 
* = R[(DrotA /D)lA(lA + 1) + (DrotB /D)lB(lB + 1)]1/2, where DrotA and DrotB are the rotational diffusion constants. â0 = a0/(4 x 82 x 82) represents the fraction of angular orientational space over which the reaction can occur, and the geometric rate is thus given by kon = 4DR x â0.
Transformation of the diffusion in a potential problem into a free diffusion problem
A crucial point in the application of equation 2 is the estimation of the angular constraints 0A, 0B, 0, and 0. We would like to estimate the ranges in mutual orientation of the two proteins for which short-range attractive forces between the atoms are sufficiently dominant to guide the two molecules into the final bound configuration, and then translate the problem of diffusional association in the attractive potential into free diffusion with an absorbing region in configurational space. To motivate this mapping, we first study two simple toy models for translational and rotational diffusion, respectively. We then use these ideas to explicitly obtain the angular constraints for real protein–protein complexes.
Toy model for translational diffusion
The reaction rate for diffusion-controlled bimolecular association of uniformly reactive spheres in the presence of a potential U(r) can be calculated exactly and is given by the expression
 | (3) |
where 
1/kT, and R1 is interpreted as the center-to-center distance between the associating partners at which the reaction is assumed to occur. Equation 3 is the classical result derived by Debye (1942). In the absence of any potential (U 0), this simplifies to the Smulochowski rate constant for free diffusion with an absorbing region of width R0, given by k(0)on = 4DR0 (in the following we use the label "0" to refer to the free diffusion problem, and the label "1" to refer to the diffusion in a potential problem).
It is clear that for R0 = R1, k(1)on > k(0)on because the presence of the (attractive) potential will increase the association rate. To find the free diffusion analog of diffusional association in an attractive potential, we increase R0 until k(0)on > k(1)on. In other words, for a given potential, we can determine the size of the absorbing region (the "capture radius") required in the case of free diffusion to obtain an association rate equivalent to that of diffusion in the potential. A similar redefinition of the effective absorbing radius to account for the presence of the potential was first introduced by Debye (1942).
A concrete illustration of this procedure is described in Figure 3 for a Lennard–Jones potential. For a broad range of parameter values, we find that it is sufficient to drop down by an energy amount of only 
(
E) = kT to enter the capture zone, corresponding to <5% of the total depth of the potential. The capture radius R0 is found to be relatively insensitive to the depth 
of the potential well, whereas its dependence on the width 
is much stronger—as it must be, because R0 is an indirect measure of the range of the potential.
|
R for the center-to-center distance of the two proteins required for the reaction to occur (see our reaction condition, equation 1), rather than using a range of allowed values (such as demanding that r 
RA + RB + R [in equation 1]), is not unreasonable.
Toy model for rotational diffusion
As another illustration of our mapping procedure, we consider two-dimensional (2D) rotational diffusion on a spherical surface in an attractive Gaussian potential U() = – exp[–( )2], with > 0. Again, we would like to translate this problem into that of free rotational diffusion with an absorbing region at = 0.
Solving the rotational diffusion equation in the presence of a potential U() yields in the diffusion-controlled limit
 | (4) |
where we let 1 
0 for U 
0 (diffusion in potential). In the case of free diffusion (U 0), equation 4 becomes:
 | (5) |
where we now choose 0 > 0. As before, we equate the association rates, equations 4 and 5, to obtain an estimate for the width 0 of the absorbing region.
For a range of parameter values, we again find that an energy drop of E U() – U(0) of the order of kT suffices to enter the capture zone (Fig. 4). We find that 0 is insensitive to the choice of but increases as expected with increasing : The range of the absorbing region reflects the range of the potential.
|
Because our toy models have only a one-dimensional (1D) reaction condition—that is, a constraint on a single degree of freedom—the question arises to what extent the relative influence of the potential on the reaction rate would change in the case of higher-dimensional reaction conditions (as used in our subsequent treatment of protein–protein interactions where we impose constraints on r, A,B, , and ). The results obtained by Zhou (1997) suggest that the influence of the interaction potential on the association rate constant is more significant for the case of two diffusing spheres bearing a circular reactive patch on each surface (i.e., where a 2D reaction condition is used for both spheres) than for the situation in which one of the spheres is taken to be uniformly reactive (i.e., where a 2D reaction condition is imposed on one sphere, but only a 1D reaction condition is used for the second sphere). Generalizing these findings, we may anticipate that an attractive interaction potential will affect reaction rates to a larger extent when the number of constrained variables in the reaction condition is increased. Because our toy models have shown that it suffices to enter the potential well by a relatively small amount to be "captured," we can conclude that for the case of a higher-dimensional reaction condition as considered in the following, an even smaller energy drop will be sufficient to enter the capture zone. From the point of view of transition state theory, our approach corresponds to identifying the transition region and then computing the flux into this region.
Mapping the protein–protein interaction funnel from the structure of a protein–protein complex
Now we apply the idea outlined above to an estimate of the angular constraints 0A, 0B, 0, and 0, needed for the application of our expression for the rate constant, equation 2. For this purpose, we have directly taken the 3D structures of the considered complexes from the Protein Data Bank (PDB).
First, the side chains of the native complexed structure were repacked by minimizing a full-atom energy function 
dominated by Lennard–Jones interactions, an orientation-dependent hydrogen-bond potential, and an implicit solvation model (Gray et al. 2003). As with all current potential functions for macromolecules, there are likely to be considerable inaccuracies in this model, but it should be emphasized that the angular constraints and rates computed here are relatively insensitive to the details of the interactions—the toy examples clearly demonstrate that once the binding funnel has been entered (which has been found to require only a small drop down in energy), the detailed form of the interaction potential has only little influence.
Second, a set of 1000 alternative structures was generated from the native complex by performing random small perturbative movements around the native conformation, and the interaction energy of these structures was evaluated using the same energy function as used in the repacking procedure. The energy landscapes defined by these alternative structures exhibit clear funnels around the native minimum.
The toy model calculations show that diffusion in such landscapes can be modeled as free diffusion with an effective "capture" region several kT into the energy funnels. To define the capture energy cutoff c below which the partners are committed to bind, we compute the average av of the energies of the five lowest-lying structures >10 Å root mean square deviation (rmsd) from the native complex (and hence outside of the native energy funnel). Because the energy cutoff cannot be determined exactly, we obtain two different estimates of the association rate setting c to either av or av – 5kT. We selected the 10 structures with the largest values of 0A + 0B in the set of structures with < c and took the averages of their values of 0A, 0B, 0, and 0 to obtain estimates for the angular tolerances used in computing the association rates.
An example is shown in Figure 5. We see that the funnel-like dependence of the energy on the rmsd and on the angular deviations is akin to the shape of the attractive potentials used in the toy models for translational and rotational diffusion discussed above. Furthermore, the location of the angular constraints resembles the position of the absorbing regions of the toy models. These similarities support our approach of deriving the angular constraints 0A, 0B, 0, and 0 from the interaction energy of perturbed protein complex structures in general, and from our method of choosing suitable energy cutoffs in particular.
|
) to compute protein–protein association rates. The effective radii RA and RB of the spheres representing the proteins are taken to be equal to the radius of gyration, Rg = (1/N)
idi (where N is the number of atoms in the protein and di is the distance of the i-th atom from the geometric center of the protein), multiplied by a correction factor of (5/3)1/2 to obtain the desired result Rg = Rs for the limiting case of a homogeneous sphere of radius Rs. The sum RA + RB is used as the value for the distance R between the centers of the two proteins at which reaction is assumed to occur; compare with equation 1. The study of our toy model for translational diffusion in the second section above has shown that this serves as a good estimate for R because the presence of short-range attractive interactions increases the effective reaction radius only slightly. The translational and rotational diffusion constants D = DAtrans + DBtrans and DA,Brot, respectively, are computed from the Stokes–Einstein relations DA,Btrans = kBT/6
RA,B and DA,Brot = kBT/8RA,B3, with = 8.9 x 10–4 Nsec/m2 (water) and T = 300 K.
Table 1 lists the 15 investigated protein–protein interactions, together with the estimated effective radii RA and RB, the angular orientational constraints 0A, 0B , 0, and 0, and the association rate constants kon determined from our theoretical expression, equation 2. For comparison, we also state the association rates obtained from a purely probabilistic model (geometric rates).
|
c = av and c = av – 5kT. For the protein complexes under study, the rates computed from the two energy cutoffs vary in average by a factor of 3, and no rates differ by more than a factor of 5 for a given complex, thus indicating the robustness of our method of estimating these rates. We observe that the angular constraints vary significantly among the investigated complexes, which suggests that our procedure of estimating these tolerances yields, indeed, characteristic and distinguishable values.
The association rate constants obtained using these angular constraints range from 104 to 106 M–1 sec–1 and are significantly higher than the corresponding geometric rates. Whereas the experimentally determined association rates of many protein–protein complexes are in this range, considerably faster rates are also observed, likely because of significant long-range interactions neglected in our model.
|
Discussion |
|---|
|
TOPAbstractIntroductionResultsDiscussionReferences |
|---|
The computed rates all lie within 104–106 M–1 sec–1, which can thus be taken as the typical diffusion-limited protein–protein association rate in the absence of attractive interactions, in good agreement with what is experimentally known for such interactions. This is several orders of magnitude higher than the geometric rate that had previously been used by various authors. Our result, therefore, shows that typical diffusion-limited association rates of proteins where no or only weak long-range interactions are present can essentially be explained with a model that is solely based on translational and rotational diffusion. Experimentally observed significantly higher rates typically suggest the presence of electrostatic steering forces, whereas much lower rates may indicate a reaction that is opposed by free energy barriers and is thus not fully diffusion-limited.
The advantage of our method over the traditional approach of BD simulations lies in the fact that our technique provides a more physically transparent insight into the resulting association rates. The differences in rates among protein complexes can be directly traced back to the sizes and shapes of the respective reactive zones in configurational space, which are determined by mapping out the binding funnel in the interaction energy landscape.
The model completely neglects possible free energy barriers caused by desolvation and/or side-chain freezing during complex formation as well as a possible slowing down of diffusion within the binding funnel caused by increased ruggedness of the landscape. Our finding that the association rates obtained with the simple diffusional model are in the range of those of many protein–protein complexes (105–106 M–1 sec–1) suggests that free energy barriers and landscape ruggedness do not have a significant impact on the dynamics of protein–protein association.
Our model provides a zeroth-order estimate of protein–protein association rates in the absence of long-range interactions. This contrasts with most previous work, which has sought to account for changes in association rates accompanying sequence changes, rather than the absolute association rate. By incorporating long-range electrostatic interactions into our diffusional model, it should be possible to develop a complete theory of association kinetics that can account for both the sequence dependence and the absolute magnitude of protein–protein association rates.
|
Acknowledgments |
|---|
The publication costs of this article were defrayed in part by payment of page charges. This article must therefore be hereby marked "advertisement" in accordance with 18 USC section 1734 solely to indicate this fact.
|
References |
|---|
|
TOPAbstractIntroductionResultsDiscussionReferences |
|---|
Debye, P. 1942. Reaction rate in ionic solutions. Trans. Electrochem. Soc. 82: 265–272.
Ermak, D.L. and McCammon, J.A. 1978. Brownian dynamics with hydrodynamic interactions. J. Chem. Phys. 69: 1352–1360.[CrossRef]
Gabdoulline, R.R. and Wade, R.C. 1997. Simulation of the diffusional association of barnase and barstar. Biophys. J. 72: 1917–1929.
———. 2001. Protein–protein association: Investigation of factors influencing association rates by Brownian dynamics simulations. J. Mol. Biol. 306: 1139–1155.[CrossRef][Medline]
———. 2002. Biomolecular diffusional association. Curr. Opin. Struct. Biol. 12: 204–213.[CrossRef][Medline]
Gray, J.J., Moughon, S., Wang, C., Schueler-Furman, O., Kuhlman, B., Rohl, C.A., and Baker, D. 2003. Protein–protein docking with simultaneous optimization of rigid-body displacement and side-chain conformations. J. Mol. Biol. 331: 281–299.[CrossRef][Medline]
Janin, J. 1997. The kinetics of protein–protein recognition. Proteins 28: 153–161.[CrossRef][Medline]
Northrup, S.H. and Erickson, H.P. 1992. Kinetics of protein–protein association explained by Brownian dynamics computer simulation. Proc. Natl. Acad. Sci. 89: 3338–3342.
Schlosshauer, M. and Baker, D. 2002. A general expression for bimolecular association rates with orientational constraints. J. Phys. Chem. B 106: 12079–12083.[CrossRef]
Schreiber, G. and Fersht, A.R. 1996. Rapid, electrostatically assisted association of proteins. Nat. Struct. Biol. 3: 427–431.[CrossRef][Medline]
Selzer, T. and Schreiber, G. 1999. Predicting the rate enhancement of protein complex formation from the electrostatic energy of interaction. J. Mol. Biol. 287: 409–419.[CrossRef][Medline]
Shoup, D., Lipari, G., and Szabo, A. 1981. Diffusion-controlled bimolecular reaction rates. Biophys. J. 36: 697–714.
Smoluchowski, M.V. 1917. Versuch einer mathematischen Theorie der Koagulationskinetik kolloider Lösungen. Z. Phys. Chem. 92: 129–168.
Vijayakumar, M., Wong, K.-Y., Schreiber G., Fersht, A.R., Szabo, A., and Zhou, H.-X. 1998. Electrostatic enhancement of diffusion-controlled protein–protein association: Comparison of theory and experiment on barnase and barstar. J. Mol. Biol. 278: 1015–1024.[CrossRef][Medline]
Zhou, H.X. 1997. Enhancement of protein–protein association rate by interaction potential: Accuracy of prediction based on local Boltzmann factor. Biophys. J. 73: 2441–2445.
CiteULike 
Connotea 
Del.icio.us 
Digg 
Reddit 
Technorati What's this?
This article has been cited by other articles:
|
|
 |
|
  M. Parat, N. McNicoll, B. Wilkes, A. Fournier, and A. De Lean Role of Extracellular Domain Dimerization in Agonist-Induced Activation of Natriuretic Peptide Receptor A Mol. Pharmacol., February 1, 2008; 73(2): 431 - 440. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 R. Alsallaq and H.-X. Zhou Energy Landscape and Transition State of Protein-Protein Association Biophys. J., March 1, 2007; 92(5): 1486 - 1502. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 N. J. Burgoyne and R. M. Jackson Predicting protein interaction sites: binding hot-spots in protein-protein and protein-ligand interfaces Bioinformatics, June 1, 2006; 22(11): 1335 - 1342. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 M. Castro, V. O. Nikolaev, D. Palm, M. J. Lohse, and J.-P. Vilardaga Turn-on switch in parathyroid hormone receptor by a two-step parathyroid hormone binding mechanism PNAS, November 1, 2005; 102(44): 16084 - 16089. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 O. Schueler-Furman, C. Wang, P. Bradley, K. Misura, and D. Baker Progress in Modeling of Protein Structures and Interactions Science, October 28, 2005; 310(5748): 638 - 642. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
0.3 kcal/mole (i.e., of order kT); this is quite insensitive to large changes in 
, a reasonable assumption for a short-range potential. Equating the resulting association rate (equation 4) with the rate for free diffusion in the presence of an absorbing region at 
