INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Diffusion of the reactants is frequently the rate-limiting step in the association of two proteins. The maximum rate constant for this process, 10⁹-10¹⁰ M¹s¹, is given by the Smoluchowski equation that describes the collision rate for two uniformly reactive spherical molecules in solution (Smoluchowski, 1917). A successful reaction between two proteins must meet the additional constraint that small reactive patches on a particular face of each protein are properly aligned. The probability of satisfying this condition in a random collision is very small, suggesting that the reaction rates should be several orders of magnitude lower than the diffusion-limited collision rates. Nevertheless, it is not uncommon to find reaction rates as high or higher than 10⁹ M¹s¹. Although this at first seems puzzling, analysis indicates that long-range electrostatic effects can heavily bias the approach of the molecules to favor reactive conditions. This effect was shown to be important for many association processes, including those of proteins with DNA (von Hippel and Berg, 1986), proteins with highly charged small molecules (Sharp et al., 1987), and proteins with oppositely charged protein substrates (Stone et al., 1989; Eltis et al., 1991; Schreiber and Ferscht, 1996; Gabdoulline and Wade, 1997; Vijayakumar et al., 1998). These systems have been thoroughly studied, and are frequently regarded as typical examples of binding phenomena.

Electrostatics is clearly not the only force that can affect the association rate. In addition to electrostatics, the most important process contributing to the binding free energy is desolvation, i.e., the removal of solvent both from nonpolar (hydrophobic) and polar atoms (Chothia and Janin, 1975). It is generally accepted that partial desolvation is always a significant contribution to the free energy in protein-protein association, and it becomes dominant for complexes in which the long-range electrostatic interactions are weak (Camacho et al., 1999). In this paper, we perform Brownian dynamics simulations to study the effects of desolvation on the rates of diffusion-limited protein-protein association.

Brownian dynamics treats each protein as a rigid body, generally a sphere, and the solvent as a viscous Newtonian liquid (Ermak and McCammon, 1978; DeLisi, 1980; Northrup et al., 1984; Northrup and Erickson, 1992; Luty et al., 1993). The method has led to a number of important results. In particular, Gabdoulline and Wade (1997) have modeled the association of the barnase-barstar complex, and found that long-range electrostatic forces alone can reconcile the high rates (10⁹ M¹s¹) observed by Schreiber and Ferscht (1996). Short-range interactions have been considered by Northrup and Erickson (1992), who have shown that a short-range locking potential can increase the association rate from 1 × 10⁵ M¹s¹ to 2 × 10⁶ M¹s¹, but did not attempt to explain the physical origin of this potential.

To establish the expected magnitude of desolvation effects on the association rates, we first perform simulations in which the interacting proteins are described by a simple model that assumes a hydrophobic interaction uniformly distributed over the entire protein surface. Results show that such short-range nonspecific interactions can significantly enhance the diffusion entrapment and thus increase the association rate. However, they also yield lengthy collisions and large nonspecific affinities that are rarely seen in real proteins (Northrup and Erickson, 1992). Therefore, we continue the analysis using a more realistic interaction potential that includes both atomic-level desolvation and short-range electrostatics.

The central assumption of the paper is that desolvation is a major contributor to the binding free energy of proteins. The magnitude of desolvation free energy is relatively well established on the basis of the free energy of transferring small molecules from water into organic solvents (Eisenberg and McLachlan, 1986, Vajda et al., 1994). We model this contribution using a structure-based atomic contact potential that has been independently validated by comparing it to various thermodynamic data, and has been shown to provide values of the desolvation free energy of proteins with remarkable accuracy (Zhang et al., 1997). Thus, there is little doubt that the desolvation force is real and, for the first time, it is included in a Brownian dynamic simulation.

Simulations are performed for two complexes in which turkey ovomucoid third domain (OMTKY) binds to -chymotrypsin and human leukocyte elastase. OMTKY has zero net charge, and, although it has higher order electrostatic multipoles, both systems represent the class of complexes in which the long-range electrostatic interactions are weak, but, under normal conditions, the binding process is still diffusion limited (Camacho et al., 1999). Typical binding rates for this type of complexes are on the order of 10⁵-10⁷ M¹s¹, about 100-fold lower than for complexes that are electrostatically steered toward the binding pocket. The simulations identify weakly specific pathways leading to a low free energy attractor of well-oriented encounter complexes, embedded in an otherwise repulsive environment. Due to these attractors, the desolvation can increase the association rate by several orders of magnitude.

As it is always the case in Brownian dynamics, the calculated absolute rates also depend on the reaction condition (Gabdoulline and Wade, 1997), and this may reduce the value of the method as a predictive tool. In our study, the simulations provide important information on the reaction condition itself. It is shown that, accounting for desolvation, the calculated and observed association rates agree if and only if the reaction condition is defined as a small ensemble of diffusion-accessible encounter complexes, suggesting that the final bound conformation is preceded by an almost unique precursor state. This interesting observation will be discussed further in the paper.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Reference coordinates and reaction condition

The receptor and the ligand are treated as spheres of radii R_r and R_l diffusing in a viscous liquid. The coordinate systems XYZ and xyz are fixed to the receptor and ligand, respectively, and their origins coincide with the corresponding centers of mass (see Fig. 1 A). The Euler angles  and  define the vector pointing to the center of the ligand, i.e., to the origin of the coordinate system xyz. Three further Euler angles _l, _l, and _l determine the relative orientation of the ligand axes xyz in the reference coordinate system XYZ.

View larger version (39K): [in this window] [in a new window]   FIGURE 1   (A) XYZ and xyz are coordinate systems fixed to the receptor and ligand centers. The position of the origin of xyz with respect to XYZ is defined by the angles (, ) and the center-to-center distance r. (l, l, l) are the Euler angles describing the ligand orientation. The optimal relative orientation of receptor and ligand is given by the angles (0, 0) and by the coordinate system x0y0z0. Notice that (0, 0) defines a vector r0 pointing to the center of the ligand, and the coordinate axes x0, y0, and z0 define the orientation of the ligand. The reaction condition is defined as a patch around this optimal point where the solid angles , x, y, and z measure the deviations from the angle (0, 0) (i.e., from the vector r0, and from the coordinate axes x0, y0, and z0, respectively). (B) The space-filling models show the van der Waals surfaces of the receptor (-chymotrypsin) and that of the ligand (OMTKY) to illustrate that these surfaces can be relatively well approximated by spheres. The overall shape and size of the second receptor considered in this work (human leukocyte elastase) is very similar to that of -chymotrypsin, and hence both molecules will be represented by spheres of effective radius Rr = 23.24 Å. The effective radius of turkey ovomucoid third domain is Rl = 14.15 Å.

The calculation of the association rate requires a reaction condition that defines the last stage of the diffusion process, after which binding would be expected to be certain. The reaction condition we use is given in terms of a reactive patch around an optimal state on each protein surface, determined by a particular position and orientation of the two molecules in contact. This orientation is given by five angles ₀, ₀; _l0, _l0, and _l0, where the first two angles determine the position of the ligand's center of mass, and the last three specify the optimal relative orientation of the ligand denoted by the axes x₀, y₀, and z₀ in Fig. 1 A. According to this condition, a simulation trajectory leads to receptor-ligand association if the receptor and ligand collide, and, at the time of the collision, the two molecules are oriented in such a way that the solid angle  around the direction (₀, ₀) in Fig. 1 A is less than a threshold _r, and, similarly, each axis of the coordinate system xyz is rotated by less than a threshold _l from the optimal axes x₀, y₀, and z₀. The relationship of this surface patch to the various reaction criteria used by other groups in Brownian dynamics simulations will be discussed further in the paper.

Brownian dynamics

The transport properties of proteins are calculated by assuming a spherical shape. Although structural asymmetries are known to yield anisotropic diffusion, it is generally accepted that, for globular proteins, these corrections should be small, and hence the spherical approximation and the resulting isotropic diffusion constants are appropriate. The translational and rotational diffusion coefficients for a sphere of radius R are given by the Stokes-Einstein relations D^trans = k_BT/6R and D^rot = k_BT/8R³, where k_B is Boltzmann constant, T is temperature and  is solvent viscosity.

According to the Brownian dynamics algorithm developed by Ermak and McCammon (1978), the time evolution of the relative displacement, r, between the reactants centers of mass of the reacting molecules is given by

(1)

where t is the time step, F is the interparticle force, and D = D_r^trans + D_l^trans denotes the sum of the unimolecular diffusion constants. S is a stochastic component of the displacement arising from random collisions of particles with solvent molecules, and is generated by taking normally distributed random numbers obeying the relationship S_k² = 2Dt. A similar expression governs the independent rotational Brownian motion. For each molecule , ( = r, l), the angular change around each of the orthogonal axes is given by

(2)

where K is the total torque on protein , and W is the stochastic term such that (W_k)² = 2D^rott. The time step t decreases monotonically from 160 ps at 500-Å separation to 0.5 ps within the desolvation layer. If a given time step leads to a protein overlap, instead of voiding the move altogether (as in Gabdoulline and Wade, 1997), we rescale t to avoid the overlap.

Throughout this paper, we study the binding of turkey ovomucoid third domain to -chymotrypsin or human leukocyte elastase. The effective radius of these molecules shown in Fig. 1 B is estimated based on the solvent-accessible surface area (Lee and Richards, 1971) of the free molecules, calculated with a water radius of 0.6 Å. The radii of chymotrypsin and leukocyte elastase differ by less than 2%, thus, for simplicity, we assume the same radius R_r = 23.24 Å for both molecules. The effective radius of turkey ovomucoid third domain is R_l = 14.15 Å. The sum R_r + R_l also agrees with the average center-to-center distance of the diffusion-accessible encounter complexes. At the temperature T = 25°C and viscosity coefficient  = 1 cP, the diffusion constants are D_r^trans = 0.0094 Å²/ps, D_l^trans = 0.0154 Å²/ps, D_r^rot = 0.0000131 rad²/ps, and D_l^rot = 0.0000578 rad²/ps.

Determination of the association rate

Brownian dynamics trajectories start with the ligand randomly placed and oriented on the surface of a sphere of radius R₀ = 60 Å in the XYZ coordinate system. The 60-Å threshold has been chosen because it is larger than the range of any intermolecular interaction. Independent trajectories are run until they either meet the reaction condition, or leave the finite diffusion space defined by a larger concentric sphere of radius R = 500 Å. The corresponding association rate constant is obtained from the equation (Northrup et al., 1984)

(3)

where  is the fraction of trajectories that satisfy the reaction condition. We note that the denominator corrects for those trajectories that reenter the diffusion space at R₀ having left at R.

Forces and torques

Let U denote the potential that determines the receptor-ligand interactions in a particular model of the receptor-ligand system. The potential depends on distance r and five angles , , _l, _l, and _l, i.e., U = U(r; , ; _l, _l, _l). Forces and torques are defined as translational and rotational derivatives of U. Using finite differences, the force component along an axis, say, X, is calculated by

(4)

where a = 0.05 Å. The torque on the ligand K_X^l is a function of the same variables, and is calculated by

(5)

where _X denotes the operator that rotates the ligand axes by  = 0.02 rad around its center of mass along the vector X. The torque on the receptor is computed by using conservation of angular momentum.

Interaction potential

We will use the binding free energy of the receptor-ligand system as the interaction potential U  G. Let G₀ denote the free energy in the unbound state. The free energy difference G = G  G₀ is calculated by

(6)

where E_coul and G_des denote the direct electrostatic (Coulombic) and desolvation contributions, respectively. In this work, we use two different desolvation models.

Model I: Uniform desolvation

Model II: Atomic-level interaction potential

Spherical approximation

Free energy landscapes

View larger version (23K): [in this window] [in a new window]   FIGURE 2   Profiles of the diffusion-accessible interaction potential of -chymotrypsin (top) and human leukocyte elastase (bottom) with OMTKY as a function of . The profiles include the low free energy attractor for both complexes. The solid lines show the potential calculated at every 1°, whereas the dashed lines correspond to the potential as used in the Brownian dynamics simulations, i.e., calculated on a coarse grid and then extended by interpolation onto the entire surface. Top:  = 85° and (l, l, l) = (20°, 180°, 180°). In the reaction condition, the reactive patch is centered at  = 85°. The arrow indicates a nearby minimum for which we also compute the association rate. Bottom:  = 80° and (l, l, l) = (20°, 300°, 40°). The reaction condition is centered at  = 95°.

Crystal structures

Three-dimensional mapping

View larger version (37K): [in this window] [in a new window]   FIGURE 3   Direct electrostatic energy (i.e., electrostatic interaction energy without including the desolvation of polar atoms) and desolvation free energy as functions of the surface-to-surface separation for four encounter complexes. Top plots are for -chymotrypsin and OMTKY; bottom plots are for human leukocyte elastase and OMTKY. Symbols × and  correspond to encounter pairs with the largest and lowest electrostatic energy, respectively. Note that, for both systems, some encounter complexes have highly unfavorable desolvation energies. Solid lines correspond to the linear approximation assumed in the Brownian dynamic simulations. Dashed lines indicate the distances defined by one, one and a half, and two layers, respectively, using a water radius of 1.4 Å.

(7)

Why do we need a spherical approximation?

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Model I: Uniform desolvation

Potential

Reaction condition

Association rates

1

1

View larger version (17K): [in this window] [in a new window]   FIGURE 4   Model I: Uniform potential. Association rate as a function of atomic solvation parameter . The number of independent runs performed to estimate the rates range from 182,000 for  = 1 to 722 for  = 10, the corresponding errors vary between 30% and 10%, respectively. Random collision rate constant is indicated by  (the solid line is the theoretical value). The + symbols denote the collision rate defined by a reactive surface patch with solid angles r = 8° around (0, 0), and 1 = 20° around x0y0z0. The * symbols show the typical time that productive pathways spend within the desolvation layer of 4.2 Å, (time scale is on the right axis).

(8)

1

1

1

1

1

1

(9)

(10)

Model II: Atomic-level interaction potential

Potential

Reaction condition

Association rates

View larger version (22K): [in this window] [in a new window]   FIGURE 5   Model II: Atomic-level interaction potential. Desolvation mediated binding: association rate for -chymotrypsin (+) and human leukocyte elastase () with OMTKY as a function of  = r = 1 (see text). A total of 40,000 and 100,000 runs were made for chymotrypsin and elastase, respectively. The errors of the calculations range from 15-30% at  = 2° to 3% at  = 90° for both chymotrypsin and elastase. The rate calculated for the same reaction condition but without any force field is indicated by the × symbols (dotted line correspond to Eq. 9). Filled squares indicate the experimental values. Inset: For comparison, we show the association rate for the lowest free-energy encounter complex in Table 1 (*), and for the local minimum indicated by an arrow in Fig. 2 () as a function of .

1

1

1

1

1

1

Association rates at nearby sites

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Background: General mechanism of association

To discuss the results of the simulations, it is useful to recall that, in the most general case, the process of diffusion-limited protein-protein association can be described by a three-step reaction mechanism,

(11)

where R denotes the receptor, L the ligand, R ··· L the nonspecific encounter pairs, and R  L the precursor state(s) leading to the docked conformation (DeLisi and Wiegel, 1981).

If long-range interactions can be neglected, the first reaction step is the random collision of the two proteins R and L, resulting in a nonspecific encounter complex R ··· L within the desolvation layer. As already mentioned, to a good approximation, the limiting rate k⁺ of this first regime is given by the Smoluchowski limit k_coll. Indeed, the overall repulsion of the force fields has little effect on k⁺. We have found that the typical lifetime of a nonspecific encounter complex R ··· L diffusing within the desolvation layer is about 4 ± 1 ns. This value is consistent with the nonspecific affinity between proteins that is estimated to be 10² M¹ or less (Northrup and Erickson, 1992).

The third reaction step in Eq. 11, i.e., the late transition between the favorable intermediate(s) R  L and the bound state RL, substantially differs from the first two steps. The onset of the late transition coincides with the need to remove steric clashes and charge overlaps in the binding mechanism. Although the first two steps are governed by diffusion, the third is a process of induced fit that requires structural rearrangements involving mostly side chains. Although the actual modeling of this regime is beyond the scope of the present paper, it is safe to assume that this late transition is not diffusive. For proteins that bind in a diffusion-controlled (or diffusion limited) reaction, the rate-limiting step must be the diffusive search for the partially desolvated intermediate(s) or precursor state(s) rather than the third step, and thus k_RL⁺ k₁.

In this paper, we focus on the kinetics of the second reaction step and on the nature of the precursor state(s) R  L. This step consists of a two-dimensional diffusive transition, driven by desolvation and short-range electrostatic forces, from a nonspecific encounter complex to the precursor state. Depending on the interaction potential, this step may or may not be important. In the well-studied association of barnase and barstar, Step 1 of the mechanism is affected by long-range electrostatic steering toward the binding site, and the encounter complex R ··· L is likely to be close to the precursor state R  L. Therefore, the role of any two-dimensional search is very limited. Without long-range electrostatics, however, any observed rate enhancement is due to short range forces that increase the probability of transition between R ··· L and R  L states.

Reaction criteria

An uncertainty inherent to any Brownian dynamics simulation is the somewhat arbitrary definition of the reaction condition. In fact, the methods described in the literature show more variation in this condition than in all the other aspects of the algorithm. A simple condition, most frequently used for model proteins, assumes that each sphere bears a small reactive patch around some point on each protein surface. The relative positions and orientations of the receptor and the ligand are defined by the direction and length of the center-to-center vector r and by the torsion angle  of rotation around r. The size of the patch on each protein is given by a solid angle  around r.

Berg (1985) and Zhou (1997) have used the above condition with no constraint on the torsion angle . It was shown that, without an interaction potential, the solid angle  = 5° yields a rate constant k_on = 6 × 10⁵ M¹s¹ (Zhou 1997). The remaining degree of freedom associated with the free rotation  has been removed by Janin (1997). If we again assume no interaction potential, we find k_on = 6 × 10⁵ M¹s¹ for  =  = 19°, whereas k_on = 1 × 10⁵ M¹s¹ requires about  =  = 14° (Janin, 1997). All these calculations account for the repeated microcollisions during the association of encounter complexes.

A contact-based reaction condition has been defined by Northrup and Erickson (1992) for spherical protein models. Each sphere with a radius of 18 Å had a set of four contact points mounted in a 17 × 17 Å square arrangement tangent to the surface. Each point on one molecule had a partner on the other molecule, and a reaction was considered to occur when N (N = 1, 2, or 3) of the four points were simultaneously within 2 Å of their partners. Notice that the cases N = 1 and N = 2 leave freedom in some orientational degrees, and hence we restrict consideration to the N = 3 case, which, apart from the 2-Å mismatch, fully specifies the relative orientations of the two molecules. Simple geometric arguments show that this mismatch can be described in terms of the surface patch model if the solid angle  and the threshold on the rotational angle  both are around 10°. Under this condition, Northrup and Erickson obtained the rate constant of k_on = 1 × 10⁵ M¹s¹ for purely diffusive association, in relatively good agreement with the results of Janin (1997). They have also shown that a short-range locking potential increases the rate constant to k_on = 2 × 10⁶ M¹s¹.

Our reaction condition, defined in terms of the vector (, ; _l, _l, _l), fully specifies the precursor state within a solid angle  = _r = _l, and, in this sense, is similar to the condition used by Janin (1997) or to the one by Northrup and Erickson (1992) at N = 3. As shown in Fig. 5, in purely diffusive association, the condition with  = 10° yields a rate constant of k_on = 1 × 10⁵ M¹s¹. Thus, despite their formal differences, the various reaction criteria generally provide similar values for the association rate without any interactions.

The nature of the precursor state

As we already mentioned, in a complex without strong long-range electrostatic interactions, one can find association rates as high as k_on = 1 × 10⁷ M¹s¹. This rate could be easily explained by using more permissive reaction criteria. For example, according to Fig. 5, the rate k_on = 1 × 10⁷ M¹s¹ can be obtained if we neglect all interactions but assume  = 24° in the reaction condition. This explanation, however, does not take into account that the short-range interactions due to desolvation definitely exist. Indeed, the thermodynamic role of desolvation is generally accepted, and its magnitude is well known. The substantial change in the desolvation free energy, occurring mostly within the first one or two water layers around the protein, necessarily yields forces associated with desolvation.

Although the mean desolvation forces are repulsive, desolvation and short-range electrostatics can substantially increase the association rate for the two complexes studied in this paper (Fig. 5). Agreement with experimentally determined rate constants can be attained under two very different assumptions. The first is assuming a very small solid angle   2°, and a necessarily fast last step to the complex RL. The second mechanism assumes a much larger ensemble of intermediates, and a slow rate to the final bound state, i.e., k_RL⁺ ~ k₁. However, this second mechanism contradicts the diffusion-limited assumption. Thus, for a diffusion-limited process, the reaction condition must correspond to the first mechanism, which involves a small set of well-defined precursor states.

We emphasize that the above results remain valid despite the potential uncertainties in the calculated rate constants. The error bars, estimated following Gabdoulline and Wade (1997), are ~30% (see Captions for Figs. 4 and 5). Comparing our results to the experimentally determined rate constants shows that agreement can be attained only by assuming a very small solid angle  in the reaction condition. In this region, the curves are so steep that even a 100% error in the rate constants has only minor effects on the required value of , and hence the conclusion concerning the nature of the precursor state is rather robust, i.e., independent of the potential errors in the calculation.

Having an almost unique precursor state suggests a relatively narrow pathway from the very restricted set of precursor states to the also unique high-affinity complex. Indeed, this late transition to the bound conformation involves going from a partially desolvated interface of 300-500 Ų to a fully desolvated interface of 1400-1600 Ų, where shape complementarity and steric effects are most important. The trajectory between these two states includes a translation of 4-6 Å, i.e., the distance separating the encounter pairs from the complex structure. It is reasonable to assume that the binding (unbinding) pathway between these two more-or-less unique states is restricted to a narrow channel in conformational space. It is only outside this channel, at partially solvated intermediates, that one could start to envisage the possibility for multiple binding (unbinding) pathways. Thus, we conclude that the late transition is a highly collimated pathway from the small set of precursor states to the complex structure. This transition should be mostly driven by a fast downhill enthalphic reduction, consistent with a locking process on the order of milliseconds (Laskowski, private communication; Lombardi et al., 1992).

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

We have performed Brownian dynamics simulations of receptor-ligand association without strong long-range electrostatic interactions. Simulations with a uniform interaction potential have shown that short-range attractive interactions substantially increase the association rate, but also result in large nonspecific affinities not seen in real proteins (Northrup and Erickson, 1992). The analysis of a realistic interaction potential, including both a well-established atomic-level desolvation and short-range electrostatics, resolves this contradiction by showing that, on average, these forces result in repulsive interactions that prevent nonspecific association. The simulations also identify weakly specific pathways leading to a low free energy attractor embedded in the otherwise repulsive environment. Along these pathways, the association rates are enhanced by the desolvation force field, both by locally increasing the diffusion entrapment and by guiding the proteins towards a small set of well-oriented kinetic intermediates. The results of Brownian dynamics simulations show that a diffusion-limited process can be reconciled with experimentally determined association rates only by assuming that the final bound conformation is preceded by an almost unique precursor state. Although desolvation can increase the association rates by several orders of magnitude, these rates are still 100-fold smaller than the ones observed for complexes in which long-range electrostatics provide the binding specificity.