Article Information

• PDF (1060 kb)
• Full Text with Large Figures
• Cited by in Scopus (14)

PubMed

• Articles by Ramzi Alsallaq
• Articles by Huan-Xiang Zhou

Related Articles

• Static and Dynamic Errors in Particle Tracking Microrheology...

  Static and Dynamic Errors in Particle Tracking Microrheology Biophysical Journal, Volume 88, Issue 1, 1 January 2005, Pages 623-638 Thierry Savin and Patrick S. Doyle Abstract Particle tracking techniques are often used to assess the local mechanical properties of cells and biological fluids. The extracted trajectories are exploited to compute the mean-squared displacement that characterizes the dynamics of the probe particles. Limited spatial resolution and statistical uncertainty are the limiting factors that alter the accuracy of the mean-squared displacement estimation. We precisely quantified the effect of localization errors in the determination of the mean-squared displacement by separating the sources of these errors into two separate contributions. A “static error” arises in the position measurements of immobilized particles. A “dynamic error” comes from the particle motion during the finite exposure time that is required for visualization. We calculated the propagation of these errors on the mean-squared displacement. We examined the impact of our error analysis on theoretical model fluids used in biorheology. These theoretical predictions were verified for purely viscous fluids using simulations and a multiple-particle tracking technique performed with video microscopy. We showed that the static contribution can be confidently corrected in dynamics studies by using static experiments performed at a similar noise-to-signal ratio. This groundwork allowed us to achieve higher resolution in the mean-squared displacement, and thus to increase the accuracy of microrheology studies. Abstract | Full Text | PDF (597 kb)

• Coarse-Grained Molecular Simulation of Diffusion and Reactio...

  Coarse-Grained Molecular Simulation of Diffusion and Reaction Kinetics in a Crowded Virtual Cytoplasm Biophysical Journal, Volume 94, Issue 10, 15 May 2008, Pages 3748-3759 Douglas Ridgway, Gordon Broderick, Ana Lopez-Campistrous, Melania Ru’aini, Philip Winter, Matthew Hamilton, Pierre Boulanger, Andriy Kovalenko and Michael J. Ellison Abstract We present a general-purpose model for biomolecular simulations at the molecular level that incorporates stochasticity, spatial dependence, and volume exclusion, using diffusing and reacting particles with physical dimensions. To validate the model, we first established the formal relationship between the microscopic model parameters (timestep, move length, and reaction probabilities) and the macroscopic coefficients for diffusion and reaction rate. We then compared simulation results with Smoluchowski theory for diffusion-limited irreversible reactions and the best available approximation for diffusion-influenced reversible reactions. To simulate the volumetric effects of a crowded intracellular environment, we created a virtual cytoplasm composed of a heterogeneous population of particles diffusing at rates appropriate to their size. The particle-size distribution was estimated from the relative abundance, mass, and stoichiometries of protein complexes using an experimentally derived proteome catalog from K12. Simulated diffusion constants exhibited anomalous behavior as a function of time and crowding. Although significant, the volumetric impact of crowding on diffusion cannot fully account for retarded protein mobility in vivo, suggesting that other biophysical factors are at play. The simulated effect of crowding on barnase-barstar dimerization, an experimentally characterized example of a bimolecular association reaction, reveals a biphasic time course, indicating that crowding exerts different effects over different timescales. These observations illustrate that quantitative realism in biosimulation will depend to some extent on mesoscale phenomena that are not currently well understood. Abstract | Full Text | PDF (638 kb)

• The Force Exerted by the Membrane Potential during Protein I...

  The Force Exerted by the Membrane Potential during Protein Import into the Mitochondrial Matrix Biophysical Journal, Volume 86, Issue 6, 1 June 2004, Pages 3647-3652 Karim Shariff, Sandip Ghosal and Andreas Matouschek Abstract The force exerted on a targeting sequence by the electrical potential across the inner mitochondrial membrane is calculated on the basis of continuum electrostatics. The force is found to vary from 3.0 pN to 2.2 pN (per unit elementary charge) as the radius of the inner membrane pore (assumed aqueous) is varied from 6.5 to 12Å, its measured range. In the present model, the decrease in force with increasing pore width arises from the shielding effect of water. Since the pore is not very much wider than the distance between water molecules, the full shielding effect of water may not be present; the extreme case of a purely membranous pore without water gives a force of 3.2 pN per unit charge, which should represent an upper limit. When applied to mitochondrial import experiments on the protein barnase, these results imply that forces between 11±2 pN and 13.5±2.5 pN catalyze the unfolding of barnase in those experiments. A comparison of these results with unfolding forces measured using atomic force microscopy is made. Abstract | Full Text | PDF (273 kb)

• …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 5, 1486-1502, 1 March 2007

doi:10.1529/biophysj.106.096024

Biophysical Theory and Modeling

Energy Landscape and Transition State of Protein-Protein Association

Ramzi Alsallaq and Huan-Xiang Zhou^,  

Department of Physics and Institute of Molecular Biophysics and School of Computational Science, Florida State University, Tallahassee, Florida

Address correspondence to Huan-Xiang Zhou, Tel.: 850-645-1336.

Interactions between proteins play central roles in diverse biological functions such as signal transduction, immune response, motility generation, and enzyme catalysis and inhibition. The mode of action is the association and dissociation of the interacting partners. The product of association is a stereospecific protein complex. Both the stability of the protein complex (measured by the association constant K_a) and the association and dissociation rate constants (k_a and k_d) are of fundamental interest. In many ways, the association process resembles the folding of a protein 1. Both are favored by short-range specific interactions, between two subunits in association while among residues within the same polypeptide chain in folding. Both are disfavored by restrictions on internal motion, i.e., relative translation and rotation for association and large-scale variations of chain conformations for folding. Great insight to protein folding has been gained from systematic studies of toy models 2,3. Here we present a study of the association process based on two types of toy models.

Echoing the protein folding process, we have previously proposed a transition state for protein-protein association 4,5. The bound state of two proteins is characterized by specific (e.g., van der Waals, hydrophobic, and electrostatic) interactions, whereas the unbound state is characterized by translational and rotational freedom. On going from the bound to the unbound, a sharp transition in interaction energy and in configurational freedom is expected. This transition serves as the outer boundary of the bound state as well as the transition state for association. A main aim of this study is to clarify the specification of the transition state.

The measured association constants for protein complexes vary from <10³ to >10¹⁵ M⁻¹. What accounts for the >10 orders-of-magnitude difference in K_a? Valuable information is provided by the structure of a protein complex determined by x-ray crystallography or NMR spectroscopy. This structure presents a representative configuration of the bound complex. However, there is no statistical mechanical reason to expect a simple relation between K_a and the interaction energy in just one representative configuration of the bound complex. The understanding of K_a is limited by uncertainties about the extent of relative translation and rotation sampled in the bound state and how the interaction energy changes with the relative motion. In recent years significant progress has been made toward a fundamental understanding of binding affinity and kinetics 1,4,5,6,7,8,9,10,11,12, but many important questions remain unresolved. With the two types of toy models studied here, in which interactions between the subunits are fully specified, we hope to address:

When two molecules associate to form a complex, the equilibrium is measured by the association constant, K_a. We now illustrate the formulation of K_a on a number of molecular models.

The development here largely follows the book of Hill 13, but with an emphasis on protein molecules. Consider a protein molecule α with internal dynamics well separated from overall translational and rotational motion. In particular, the separability of internal and overall motion is the basis of the model-free approach to the analysis of NMR relaxation data 14. Let the overall translation be described by the displacement vector r_α, the overall rotation be described by an orientation vector ω_α, and the internal degrees of freedom be represented by a vector x_α. If the potential of mean force, after considering the solvent degrees of freedom, of the protein molecule is U_α(x_α), then the configurational integral is

(1)

Now suppose that there are three species, protein A, protein B, and their complex C. For protein A or B, Eq. (1) gives the configurational integral when the subscript α is replaced by A or B. For the complex, translations of the subunits can be recast as the overall translation R plus a relative translation r=r_B-r_A. The potential energy of the complex can be written as U_C(r, ω_A, ω_B, x_A, x_B), and its configurational integral is

(2)

The equilibrium constant for the association of proteins A and B to form the 1:1 complex C is then 13

(3a)

(3b)

(4)

(5)

So far, the individual translations of the subunits in the complex have been recast as the overall translation of the complex plus relative translation within the complex, but the individual rotations of the subunits have been retained in the formulation of K_a. One way to separate overall and relative rotations is to 1), select a fixed orientation of protein A and then sample different orientations of protein B; and 2), rotate protein A and protein B together to different orientations. Let the orientation of protein B relative to the selected orientation of A be ω, and the orientation of the subunits together be Ω, then

(6)

It has often been said that, when two proteins form a complex, six translational and rotational degrees of freedom are lost. This statement is misleading as it neglects the relative translational (r) and relative rotational (ω) motion within the complex. A completely rigid complex, i.e., one without any relative translation or rotation, has an association constant that is given by an integral over a single point, which is zero unless the interaction potential well is infinitely deep.

The formulation of the association constant presented above is based on separating the relative translational and rotational degrees of freedom (r and ω) from the internal degrees of freedom (x_A and x_B) of the subunits. Tidor and Karplus 15 used normal-mode analysis to study the contribution of relative motion within the bound complex. The advantage of their approach is that there is no need to explicitly separate r and ω from x_A and x_B. The advantage of our formulation is that there is no need to assume the potential U(r, ω_A) as harmonic.

The change in chemical potential upon the association of proteins A and B is 13

(7a)

(7b)

(7c)

Consider two spherical proteins, with an interaction energy U(r) depending on the interprotein distance r. A complex is considered formed when r is within an outer limit r^‡ defining the bound state. For this model, the association constant is given by 5,6,11

(8)

(9)

With the individual rotations cast aside, the two spherical subunits can be viewed as point masses, and the complex as a diatomic molecule. While each subunit before association has three translational degrees of freedom, the complex has three degrees of translation for the center of mass, two degrees of rotation (around two perpendicular axes through the center of mass), and one degree of vibration for the intersubunit distance. If the masses of the subunits are m_A and m_B, then their partition functions for translational motion are

(10)

(11)

In this simple model, individual rotations of the subunits and their relative translation are uncoupled, and the relative translation can be further separated into rotation of the complex and vibration within the complex. Rotation of the complex is equivalent to a change in the direction of the relative translation, whereas vibration within the complex is just a change in the magnitude of the relative translation. For more complicated models, individual rotations and relative translation become coupled. It then becomes impractical to introduce harmonic approximations to any of the translational and rotational degrees of freedom.

Equation (6) can be interpreted as the expression for the association constant of two proteins modeled as rigid bodies that are interacting with an interaction potential U(r, ω). One expects U(r, ω) to have a deep minimum, which identifies the bound state. As already alluded to, the bound state is not just a single configuration of the complex. Configurational sampling around the energy minimum must take place in the bound state, and the integral in Eq. (6) should reflect this sampling. It is clear that K_a is not solely determined by the minimum interaction energy, U_m, in the bound state. Equally important are the variation of the interaction energy with relative motion and the size of the configurational space of the bound state. Unfortunately, x-ray crystallography cannot tell the extent of configurational space sampled by the constituent proteins in the bound state. Some relative motion around the energy minimum is prevented by steric clashes between the subunits. Steric clashes thus set an inner boundary for the bound state. The outer boundary remains to be specified.

In this study we specify the six relative translational and rotational degrees of freedom of the bound state in the following way (see FIGURE 1A). Protein A is fixed in space, with the center of the binding site at the origin of a laboratory-fixed coordinate system. Protein B is allowed to translate and rotate. The position of the binding site of protein B in the laboratory-fixed coordinate system is identified as the relative translation vector r. With the magnitude of the r denoted as r and the direction of r specified by polar and azimuthal angles θ and ϕ, dr=r²drsinθdθdϕ=−r²drdcosθdθdϕ. The orientation ω of protein B in the laboratory-fixed coordinate system consists of a body-fixed unit vector e (with the polar and azimuthal angles ξ and ζ) and a rotation angle χ around the unit vector. In terms of these angles, dω=sinξdξdζdχ=−dcosξdζdχ.

Display large version of this figure

FIGURE 1

(A) Definitions of the three translational (r) and three rotational (e and χ) degrees of freedom. (B) The hemisphere and (C) crater models in their minimum-energy configurations. (D) Side view of the interface between two subunits with two types of interaction sites: filled for h and open for p.

Equation (6) can be rewritten as 4

(12)

The total sampling volume in the six-dimensional configurational space is V₀=16π³(1-cosθ₀)r₀ when r, cosθ, ϕ, cosξ, ζ, and χ are uniformly sampled. Along r, there is a geometric factor r², which should appear as a weighting factor when the uniform sampling is used to calculate averages. Specifically, r, the average of a quantity A, 〈A〉, is calculated as 〈r²A〉′/〈r²〉′, where 〈⋯〉′ means averaging with uniform sampling. Of all the configurations distributed within the sampling volume, some do not contribute to K_a because they involve steric clashes between the subunits. Let the fraction that avoids clashes be f_c. Within this fraction, only a subfraction (f_b) is in the bound state. The configurational volume of the bound state is given by

(13a)

(13b)

(14)

(15)

We implement two type of toy models to illustrate the calculation of the association constant. They have different interface shapes, mimicking in a small measure the great variety of interface shapes of actual proteins. The first, called the hemisphere model, consists of two matching hemispherical proteins (both with a radius denoted by R), which form a whole sphere in the bound state (see Fig. 1). The binding site on each subunit is a flat circle with an area S=πR². The second toy model, called the crater model, consists of a spherical protein with a crater to which another spherical protein snugly fits in the bound state (see Fig. 1). The radii of the two proteins are denoted as R_A and R_B, respectively. In this case the binding site is curved on each side. If the polar angle spanned by the binding site on protein B is γ, then the interface area is . In this study cosγ is set to 0.5, so the interface areas of the two models are identical when R=R_B. In both models, the body-fixed unit vector e on protein B is chosen to be the normal vector located at the center of the binding site, pointing toward the interior of the protein. The coordinate systems for r and ω are defined such that the configuration in which the two subunits are perfectly matched corresponds to r=0, cosθ=cosξ=1, and χ=0. Note that for the hemisphere model, when r<R, cosθ must be greater than 0 to avoid collision between the two subunits.

To model interactions, matching loci on the binding sites of the two subunits are randomly selected, with a minimum separation of s_m=3.5 Å among the loci on either side. To ensure stereospecificity of the bound complex, there are two types of interaction loci (labeled h and p). The total numbers of h and p loci are denoted n_h and n_p, respectively. Each locus on protein A potentially interacts with all the loci on protein B and vice versa. The interaction energy between two loci across the interface is a square well (between two h-loci or two p-loci) or square barrier (between an h- and a p-locus), with a width r_w=3.5 Å (see FIGURE 1D). When the two subunits collide, the interaction energy is infinite; otherwise it is given by

(16)

Two sets of parameters are implemented for both the hemisphere model and the crater model. In the first set, R=20 Å for the hemisphere model while R_A=25 Å and R_B=20 Å for the crater model. The interface areas of the two models are both 1257 Ų. Over this area 30 interaction loci are distributed, of which 18 were h-loci and 12 are p-loci. In the second set, R=17 Å for the hemisphere model, while R_A=21 Å and R_B=17 Å for the crater model. The interface areas of the two models are now 908 Ų. Over this area 20 interaction loci are distributed, of which 12 are h-loci and 8 are p-loci. The interaction locus densities are ∼1 in every 45 Ų of interface area for all models. For easy reference, the hemisphere and crater models with the larger surface area will be denoted as HL and CL, respectively, and the corresponding ones with the smaller interface area will be denoted as HS and CS, respectively. The sampling bounds r₀ and cosθ₀ for the models are listed in Table 1,Table 2.

Table 1 Transition state and bound state properties of the toy models with an interface area of 1257 Å2

	VariablesHemisphere model	Crater model
	Sampling range of coordinate
	r0 (Å)	10	6	10	6
	cosθ0	0	0	−1	−1
	Mean and standard deviation of coordinate in transition state
	(Å)	3.3±0.7	3.4±0.5	3.7±0.9	4.0±0.7
		0.8±0.2	0.8±0.2	0.7±0.2	0.7±0.2
		3.0±1.8	3.0±1.8	2.8±1.9	3.0±2.0
		0.995±0.004	0.994±0.005	0.988±0.011	0.987±0.009
		2.9±1.7	2.9±1.7	3.6±1.6	3.4±1.7
		0.004±0.15	−0.001±0.19	−0.01±0.17	0.003±0.20
	Energetic and geometric parameters of bound state
	−βUm	54.98	56.55	54.98	54.98
	−βU‡	23.56	21.99	25.14	21.99
	〈exp(−βU)〉b	6.79x1020	6.26x1020	4.44x1020	3.01x1020
	−βGint	47.97	47.89	47.54	47.15
	fc	7.18x10−3	2.52x10−3	1.66x10−3	0.510x10−3
	fb	2.26x10−3	11.96x10−3	2.10x10−3	14.57x10−3
	〈r2〉′b (Å2)	6.75	7.13	7.54	8.78
	(10−3 Å3)	6.88	8.10	3.31	4.92
	Ka (1015 M−1)	2.8	3.1	0.89	0.89

Table 2 Transition state and bound state properties of the toy models with an interface area of 908 Å2

	VariablesHemisphere model	Crater model
	Sampling ranges of coordinate
	r0 (Å)	10	6	10	6
	cosθ0	0	0	−1	−1
	Mean and standard deviation of coordinate in transition state
	(Å)	3.4±0.8	3.5±0.7	3.5±0.6	3.3±0.5
		0.7±0.2	0.7±0.2	0.7±0.2	0.7±0.2
		3.0±1.9	3.0±1.9	3.0±1.8	3.1±1.8
		0.994±0.006	0.993±0.007	0.986±0.010	0.988±0.008
		3.5±2.1	3.5±2.1	3.3±1.9	3.2±1.9
		−0.004±0.23	−0.002±0.25	−0.03±0.18	−0.03±0.14
	Energetic and geometric parameters of bound state
	−βUm	39.27	39.27	36.13	36.13
	−βU‡	18.85	17.28	17.28	18.85
	〈exp(−βU)〉b	1.36x1014	0.953x1014	0.111x1014	0.154x1014
	−βGint	32.54	32.19	30.04	30.36
	fc	10.1x10−3	3.54x10−3	2.47x10−3	0.770x10−3
	fb	2.46x10−3	14.13x10−3	2.09x10−3	9.65x10−3
	〈r2〉′b (Å2)	6.55	7.23	7.37	6.58
	(10−3 Å3)	10.2	13.6	4.78	3.69
	Ka (108 M−1)	8.4	7.8	0.32	0.34

The outer boundary of the bound state dictates the rate constant at which the complex is formed by diffusion. Experimental data on the diffusion-controlled association rate thus provide valuable information for specifying the outer boundary. In the absence of long-range interactions, the diffusion-controlled association rate constant, , is typically in the range of 10⁵–10⁶ M⁻¹ s⁻¹17. The parameters used in the models should lead to values of that are in this range. We carry out this important check by calculating for each model through Brownian dynamics simulations 18.

The outer boundary of the bound state plays a critical role in the transition-state theory for the protein-protein association 4,5,17. In this theory, the outer boundary is taken to be the transition state, and the association rate constant in the presence of long-range electrostatic interactions is calculated as

(17)

The toy models are designed to capture two essential properties of protein complexes:

Configurations of the subunits in each model system are randomly generated, with the relative separation (r) restricted to within 10 Å. Each configuration is checked for collision between the subunits. If no collision occurs, the interaction energy is calculated. Typically energy calculations are made on 10million configurations. FIGURE 2A displays a scatter plot of the interaction energy versus the rotation angle χ for the HL model. A striking feature of the plot is the sudden transition between the bound state in which relative rotation is restricted and interaction between the subunits is strong and the unbound state in which relative rotation is unrestricted but interprotein interaction is weak. This contrast is manifested by the standard deviation of χ (σ_χ) sampled at different energy levels (FIGURE 2B). The transition can be conveniently located by the parameter

(18)

Display large version of this figure

FIGURE 2

(A) Scatter plot of the interaction energy versus the rotation angle χ for the HL model. For clarity, the full range of χ is evenly divided into 500 bins and each bin contains at most one sampled χ-value at each energy level. The total number of sampled configurations is 10⁷ (r₀=10 Å). The transition-state energy level is indicated by dark points. (B) The standard deviation of χ and the parameter (Eq. (18)) at different energy levels. An arrow indicates the transition-state energy level, where is maximal.

Another feature of the U-versus-χ scatter plot is the profile of the lower bounds of the sampled energies at different χ-values. For example, the lower bounds at ±90° are much lower than those at ±120°. The lower bounds in U are reached when the other five coordinates are close to those in the perfectly matched configuration. Starting with the perfectly matched configuration (in which χ=0°), as the value of χ is changed toward ±180°, the energy shows nearly the same profile as in FIGURE 2A. This profile is a reflection of the distribution of the interaction loci within the binding sites.

The energies and statistical distributions in χ and the other coordinates for the transition states of the four toy models are listed in Table 1,Table 2. The span of allowed values of the relative separation r experiences a sharp transition, as illustrated by a scatter plot in FIGURE 3A for the HL model, similar to the situation for the rotation angle χ. The transition also occurs at the same energy level U^‡ as found in the χ-dependence. The mean values of r among the transition-state configurations of the four toy models all fall within 3 to 4 Å; the standard deviations are between 0.5 and 1 Å.

Display large version of this figure

FIGURE 3

(A) Scatter plot of the interaction energy versus the relative separation r for the HL model. The full range of r (0–10 Å) is evenly divided into 500 bins and each bin contains at most one sampled r-value at each energy level (out of a total of 10⁷ sampled configurations). (B) The free-energy functional W(r) and its energetic and entropic components.

In contrast to χ and r, sampling along the other four coordinates shows much less variation between the bound and unbound states. Avoidance of collision forces the direction of the translation of subunit B to be away from subunit A (i.e., cosθ>0). This is so even for the crater models, for which the lower bound of cosθ is set to −1. There are no significant differences in the ranges of cosθ and ϕ sampled between the bound and unbound states. Avoidance of collision also significantly restricts the direction of the body-fixed unit vector e. Only values of cosξ>0.8 are sampled in each of the four models (freedom in cosθ and cosξ will eventually be regained when r is greater than the sum of the subunit radii), with values in the bound state restricted to ∼0.98, corresponding to a polar angle of ∼10°. Because of the restriction on cosξ, there is apparent freedom in the azimuthal angle ζ of the unit vector e (note that the value of ζ is irrelevant at cosξ=1).

Though it is not possible to sample the full span of possible r values (0 to ∞), the regions of interest, i.e., the bound state and the transition to the unbound state, occur well below the upper bound r₀=10 Å and hence are well sampled. To make sure that conclusions are not influenced by the specific value of r₀, the configurational spaces of the four toy models are also sampled with r₀=6 Å. All the main results presented above are confirmed, with a possible small shift in U^‡ (Table 1,Table 2).

Along the r coordinate, the interplay between interaction energy and configurational freedom in the transition from the bound state to the unbound state can be elucidated by separating the Boltzmann weight into energetic and entropic contributions. The energy function E(r) can be defined from the average Boltzmann factor among the N(r) allowed configurations in a bin [r-Δr/2, r+Δr/2]:

(19a)

(19b)

(19c)

FIGURE 3B displays the free energy functional and its energetic and entropic components for the HL model. Note that, even at r=10 Å, βE(r) and S(r)/k_B are still quite significant (at ∼−10 and −4, respectively). The interactions contributing to E(r) at such large separations between the centers of the binding sites come from loci on the peripheries of the binding sites. These interactions may hold the subunits together to allow them time to search for the bound state. Such an “entrapment” effect has been seen in Brownian dynamics simulations 23. The energy and entropy functions have very different dependences on r. E(r) decreases sharply as the subunits enter the bound state (at around the mean r-value, 3.4 Å, of the transition-state configurations). The change in S(r) is more gradual. The asynchronous changes in E(r) and S(r) do not seem to lead to a significant free-energy barrier, unlike what was speculated previously 5. How the change in internal degrees of freedom during the association process affects the energy and entropy functionals remains to be studied.

A similar free energy functional, W(r, χ), which depends both on r and χ, can also be defined. The energy component is again given by the average Boltzmann factor according to Eq. (19a), but with N(r) replaced by N(r, χ), the number of allowed configurations within a two-dimensional grid with r in [r-Δr/2, r+Δr/2] and χ in [χ-Δχ/2, χ+Δχ/2]. The entropy component is given by Eq. (19b), but with N(r) in the numerator replaced by N(r, χ) and an additional factor, 2π/Δχ, is inserted. FIGURE 4AD, displays W(r, χ) for the four toy models. The functional presents a funnel-like energy landscape, with the deep well of the bound state surrounded by a broad shallow basin.

Display large version of this figure

FIGURE 4

The free-energy functional W(r, χ) for the (A) HL, (B) CL, (C) HS, and (D) CS models.

Our aim is to extend the study of protein-protein association from toy models to actual protein-protein complexes. In that case, calculation of interaction energies based on realistic molecular models becomes a formidable challenge. Therefore we have sought alternatives to the interaction energy for obtaining the energy landscape of protein-protein association.

Inspired by the use of contacts in studying protein folding, we have tested different contact-based choices and found a reasonable surrogate in N_c, the sum of native contacts and nonnative contacts. The former are taken as formed by cognate pairs of interaction loci within a distance of r_w=3.5 Å whereas the latter are taken as formed by noncognate pairs of interaction loci within a shorter distance threshold r_w′=2.5 Å. Scatter plots of N_c versus the rotation angle χ and the relative separation r are shown in FIGURE 5AB, for the HL model. These pictures are qualitatively very similar to those found for the interaction energy (see FIGURE 2AA and FIGURE 3AA). In the N_c-versus-r scatter plot, a void appears at the r=0 and N_c=0 corner, because when the separation is small the two subunits will inevitably make at least a few contacts. The transition-state N_c level, , can be found from the maximum of the parameter (N_c) that is defined analogous to Eq. (18). The value of thus found is 20, 21, 15, and 14, respectively, for the HL, CL, HS, and CS models.

Display large version of this figure

FIGURE 5

Scatter plots of the total contact number versus (A) the rotation angle χ and (B) relative separation r for the HL model. Dark points indicate .

The transition-state configurations obtained with N_c are very similar to those obtained with energy. The average energies of the N_c-based transition-state configurations are 24, 25, 23, and 18, respectively, in units of k_BT. These are close to the values of βU^‡ listed in Table 1,Table 2. In addition, the means and standard deviations of the six coordinates in the N_c-based transition-state configurations are close to their energy-based counterparts (data not shown).

A “free-energy functional” W_Nc(r) can be defined analogous to W(r), with −βU in Eq. (19a) replaced by N_c. Note that the entropic components of W_Nc(r) and W(r) are identical. The energetic component E_Nc(r) shows high correlation with its counterpart E(r) (e.g., with R²=0.97 for the HL model). Similarly, in analogy to W(r, χ), a two-dimensional free energy functional, W_Nc(r, χ), can be defined. FIGURE 6A displays this functional for the HL model, which presents the same funnel-like energy landscape as seen in FIGURE 4A. The correlation between E_Nc(r, χ) and E_Nc(r, χ), when the latter is <−15 k_BT, is shown in FIGURE 6B for the HL model.

Display large version of this figure

FIGURE 6

(A) The free-energy functional W_Nc(r, χ) for the HL model. (B) Correlation of the energetic component E_Nc(r, χ) with its counterpart E(r, χ). The values of the latter is restricted to <−15 k_BT. This range of values covers the bound state and much of the transition region to the unbound state (correlation outside this range deteriorates). The line of linear regression (with zero intercept) is shown.

The formalism developed for the toy models can be directly applied to actual protein-protein complexes. For illustration, here we present the application to the barnase-barstar complex, a system that has been studied extensively 4,5,19,21,24,25,26,27,28,29,30. We treat each subunit as a rigid body. The binding sites on the two subunits are identified by heavy atoms making interfacial contacts which are <5 Å in the x-ray structure of the complex (Protein Data Bank entry 1brs chains C and F; Fig. 7) 25. There are a total of 109 such atoms on the barnase side and 101 on the barstar side. The geometric center of this collection of interface atoms and the normal vector of their least-square plane are used to define the coordinate systems for r and ω. The laboratory-fixed coordinate system has its origin at the geometric center and its z axis along the normal vector. Barnase is then fixed in this coordinate system. The geometric center and the normal vector is body-fixed on barstar, which is then allowed to translate and rotate. The position of the geometric center fixed on barstar defines r; the normal vector fixed on barstar becomes the body-fixed unit vector e, which together with a rotation angle χ define ω. The x-ray structure corresponds to r=0, cosθ=cosξ=1, and χ=0.

Display large version of this figure

FIGURE 7

The interaction-locus atoms of the barnase-barstar complex. Barnase and barstar are shown in blue and gray, respectively; the interaction-locus atoms are shown as blue or red spheres.

Interaction loci are identified from the collection of interface atoms. For each interface atom, the shortest cross-interface contact with a heavy atom in the x-ray structure is found. All such cross-interface contacts are then sorted in ascending order of contact distances. If a contact-forming atom is within s_m=3.5 Å of an atom on the same protein which forms a shorter contact, the longer contact is eliminated from the list. In the end, a total of 17 distinct contacts are retained. Fig. 7 displays the interaction loci, i.e., atoms forming those cognate contacts. For each interaction-locus atom, a contact radius is defined as half of the contact distance with its partner.

For configurations that do not involve steric collision between the two subunits, the number of contacts, N_c, is found by summing the number of native and nonnative contacts. A native contact is formed by an interaction-locus atom with its cognate partner, with a distance that is not longer than the value found in the x-ray structure by r_w=3.5 Å. A nonnative contact is formed by a noncognate pair of interaction-locus atoms, with a distance that is not longer than the sum of their contact radii by r_w′=2.5 Å. In particular, in their x-ray structure, there are 21 nonnative contacts in addition to the 17 native contacts, resulting in N_c=38. Steric collision is detected whenever a pair of atoms on the two subunits is closer than a collision distance. For the purpose of detecting collision, atoms are classified into three types: hydrogen, polar (nitrogen and oxygen), and nonpolar (carbon and others). The collision distance within one type or between two types of atoms is set to the minimum distance of such contacts in the x-ray structure of the complex. The resulting collision distances are: 2.64 Å between polar atoms, 3.48 Å between nonpolar atoms, 3.11 Å between polar and nonpolar pairs, 2.14 Å between hydrogens, 1.63 Å between polar and hydrogen atoms, and 2.51 Å between nonpolar and hydrogen atoms.

Scatter plots of N_c versus the rotation angle χ and the relative separation r are shown in Fig. 8 for the barnase-barstar complex. These plots show resemblance to corresponding plots for the toy models (see Fig. 5). The transition from the bound to the unbound state occurs at N_c=14 for configurations sampled with r₀=6 Å. The mean and standard deviation of χ in the transition-state configurations are 0.01 and 0.31 radians (or 0.6° and 18°), respectively. The mean and standard deviation of r in the transition state are 4.9 and 0.5 Å, respectively. These and other statistics of the transition state are collected in Table 3. Representative configurations in the transition state are shown in Fig. 9. In Fig. 10 we display the free-energy functional W_Nc(r, χ). Similar to the situation found in the toy models (see FIGURE 6A), this functional exhibits a funnel leading to the minimum-energy configuration (as given by the x-ray structure of the complex).

Display large version of this figure

FIGURE 8

Scatter plots of the total contact number versus (A) the rotation angle χ and (B) relative separation r for the barnase-barstar complex. Dark points indicate . To ensure adequate sampling of the whole range of r from 0 to 10 Å, three independent runs are carried out with the upper bound of r set to 4, 6, and 10 Å, respectively. The sampled configurations are then combined.

Table 3 Transition state and bound state properties of barnase-barstar complex

	r0 (Å)	6
	cosθ0	−1
	(Å)	4.9±0.5
		0.89±0.08
		3.2±2.5
		0.92±0.05
		2.8±1.7
		0.01±0.31
	Nc;max	30
		14
	fc	0.101x10−3
	fb	0.059
	〈r2〉′b (Å2)	19.5
	(10−3 Å3)	8.75
	Interface area (Å2)*	797
	Ka (1012 M−1)†	8

* Taken as half of the buried solvent-accessible area calculated with a 1.4 Å probe radius. † Schreiber and Fersht 24. The value listed is for an ionic strength of 125mM.

Display large version of this figure

FIGURE 9

Representative transition-state configurations. Barnase, shown in blue, is fixed in the laboratory frame whereas barstar is allowed to translate and rotate. The body-fixed unit vector e on barstar in different configurations is shown as arrows. For one particular configuration, the arrow is in blue and the corresponding structure of barstar is shown in gray.

Display large version of this figure

FIGURE 10

The free-energy functional W_Nc(r, χ) for the barnase-barstar complex. Regions not covered by the free-energy surface are not sampled.

While the transition from the unbound to the bound state is qualitatively similar to those in the toy models, a major difference in the barnase-barstar complex is that the interface involves atomic details. Because of this, when the relative separation is small, collision between the subunits can be avoided only if they are nearly aligned for complex formation. In particular, at small r rotation around the body-fixed unit vector becomes very restricted, and only a narrow range of χ around 0° can be sampled. This explains why in Fig. 10 the W_Nc(r, χ) surface does not cover the full sampled range of χ.

Comparing FIGURE 8A with FIGURE 5A, it can be seen that fluctuations in the lower bounds of the sampled N_c values across the range of χ are much more prominent in the barnase-barstar complex. In the toy models, the variation of the lower bounds of N_c with χ reflects the distribution of the interaction loci within the interface. Unlike in the toy models, a change in χ in the barnase-barstar complex can lead to collision between the proteins. Avoidance of collision thus introduces additional variations in the lower bounds of N_c. Alignment of the proteins at small separations required by avoidance of collision also explains the expanded void at the r=0 and N_c=0 corner in the N_c-versus-r scatter plot (comparing FIGURE 8B against FIGURE 5B). As the proteins are more aligned, more contacts will also form.

As described in the previous section, the outer boundary of the bound state is specified as the transition region in translational and rotational freedom on going from the bound to the unbound state. This specification allows for an unambiguous calculation of the association constant from the integration in Eq. (6). On the other hand, the numerical value of K_a should not be sensitive to the precise specification of the outer boundary of the integration.

Table 1,Table 2 list geometric and energetic information for the bound state of the toy models. The operational definition of the bound state is given by the energetic criterion U≤U^‡. As noted in the previous section, the relative separation and orientation between the two proteins in the bound state are severely restricted. Specifically, the separation distance r is restricted to within ∼4 Å, the direction of the body-fixed unit vector e on protein B is restricted to within ∼10°, and rotation around this vector is restricted to within ∼15°. The resulting configurational volume of the bound state is ∼10⁻² ų. The free energy of interaction in the bound state, G_int, defined in terms of the average Boltzmann factor via Eq. (14), is ∼5–10 k_BT higher than the minimum-energy (U_m) found from configurational sampling. The association constant is 2.8x10¹⁵ and 0.89x10¹⁵ M⁻¹, respectively, for the HL and CL models and 8.4x10⁸ and 0.32x10⁸ M⁻¹, respectively, for the HL and CL models.

To illustrate the insensitivity of the value of K_a to the precise specification of the outer boundary of the bound state, in Fig. 11 we plot K_a for the HL model calculated at different levels of βU^‡. It can be seen that essentially the same value of K_a is obtained from the sampled configurations as long as βU^‡>−48. When the bound state is located in a deep energy well, the integral of the Boltzmann factor for evaluating K_a (Eq. (6)) is dominated by a small region around the energy minimum, and the precise specification of the limits of the integral has no consequence on the numerical value of K_a. This is the same reason why the method of steepest descent works so well for evaluating integrals of functions with sharp maxima. However, as the sampling volume increases, eventually the integration and the resulting K_a value go to infinity.

Display large version of this figure

FIGURE 11

Calculated value of K_a versus the energy level used to define the outer boundary of the bound state.

The values of K_a calculated with configurations sampled with r₀=10 and 6 Å are essentially identical. This is despite the fact that in the latter configurations the collision-free fraction f_c is lower by approximately threefold and the bound fraction f_b is higher by approximately sixfold. Moreover, the configurational volumes of the bound state obtained from the two sets of sampling agree closely. For the HL model, the minimum interaction energy obtained from the r₀=6 Å configurations is lower by one level, indicating that there is uncertainty in obtaining the absolute minimum interaction energy by sampling.

While the calculation of the association constant using different specifications of the outer boundary of the bound state confirms that the numerical value of K_a is unequivocally determined, it does indicate that the further breakup of K_a into the average Boltzmann factor and the configurational volume (Eq. (12)) is to a certain extent arbitrary. On the other hand, the diffusion-controlled rate to reach the bound state is sensitive to the size of the configurational volume, hence data on serve as a determinant of .

We have carried out Brownian dynamics simulations to calculate the diffusion-controlled rate constant to reach the bound state, using an algorithm developed previously 18. With the outer boundary of the bound state as specified in Table 1,Table 2, is found to vary from 0.6x10⁵ to 3.7x10⁵ M⁻¹ s⁻¹ for the four models. These values fall within the range of 10⁵ to 10⁶ M⁻¹ s⁻¹ observed experimentally for the diffusion-controlled association of proteins in the absence of long-range electrostatic enhancement 17. Thus our specification of the outer boundary of the bound state appears reasonable. Again, the configurational volume of the bound state thus obtained is ∼10⁻² ų.

In the previous section it is seen that, according to statistics of the six translational and rotational coordinates of sampled configurations, the transition state obtained with the number of contacts, N_c, closely mimics that obtained with energy. The configurational volume of the bound state calculated with the two approaches also agrees to within a factor of 3, as does , the diffusion-controlled rate constant for reaching the bound state.

In theoretical studies of the effects of mutations on the association constant, a common practice is to calculate mutational effects on the interaction energy in the x-ray structure of the complex 28,31,32. This is equivalent to the approximation

(20)

Ten single, double, triple, quintuple, and decuple mutations each are made by deleting 1, 2, 3, 5, or 10 interaction loci on subunit A of the CL model. For the single mutations, the mutation loci are selected randomly. For the other mutations, the most closely clustered sets of loci are selected. The association constant for each mutant is then calculated in the same way as the original (“wild-type”) CL model.

For the wild-type CL model, 38 configurations are found to have the minimum energy U_m=−54.98 k_BT. The mutations are applied to these configurations, and the changes in the interaction energy (from U_m