INTRODUCTION

Top Abstract Introduction Methods Results & Discussion References

More than half a century ago, Zimm (1946) stated that "it is a task of statistical mechanics to relate the thermodynamics of solutions to the properties of the molecules that compose them." He used his theoretical results in part to interpret thermodynamic data on protein solutions. With the advent of modern biotechnology, the incentive for such analyses of protein solutions has become much stronger, especially for application to downstream processing operations (Haynes et al., 1993; Watanabe et al., 1994; Coen et al., 1995). Protein thermodynamics in vivo is an even greater challenge, being confounded by the complexity of the intracellular environment. However, it is one that is likely to receive more attention as rapidly emerging genome sequences are used to generate protein structural information and ultimately to develop morphological and functional pictures of intact organisms (Bassingthwaighte, 1997; Wilkins et al., 1997). Models describing such systems should faithfully represent the underlying physics and not just fit experimental data. The purpose of the present work is to contribute to such an understanding.

The thermodynamic property of proteins examined by Zimm (1946) was the osmotic second virial coefficient B₂₂, defined by

(1)

where  is the osmotic pressure, c_p is the protein concentration (in mass units), R is the gas constant, T is the absolute temperature, and M_w is the protein molecular weight. B₂₂ reflects the size and direction of deviations from ideality of the solution, with an ideal solution obeying the van't Hoff relation, i.e., Eq. 1 with B₂₂ and higher order terms set to zero. At the molecular level, B₂₂ reflects the nature of protein-protein interactions: Eq. 1 shows that a positive value of B₂₂ indicates an osmotic pressure larger than that for an ideal solution, which can be interpreted intuitively to reflect predominantly repulsive interactions among protein molecules, with the converse true for negative B₂₂.

The insights that B₂₂ can provide into molecular interactions have driven ongoing experimental work and parallel efforts to model the interactions. A more complete understanding of the molecular origins of B₂₂ can shed light on experimental means of manipulating the thermodynamic properties and phase behavior of proteins and may offer new approaches for extracting partial structural information from solution property measurements. The quantitative relation of B₂₂ to molecular interactions is usually presented in terms of the orientationally averaged potential of mean force (PMF), W(r₁₂), where r₁₂ is the intermolecular center-to-center distance (McMillan and Mayer, 1945; Hill, 1987):

(2)

Here B'₂₂ = B₂₂M_w², and W is essentially the interaction free energy between two protein molecules in an averaged solvent environment, and k is the Boltzmann constant. (Some care is necessary in comparing calculated B₂₂ values with experiment because of differences in the thermodynamic conventions used (Cabezas and O'Connell, 1993), but the effects are not large (Coen et al., 1995).) The portion of the integral corresponding to overlap, where W  , leads to the excluded volume contribution to the virial coefficient, which is always positive. Remaining features of W(r₁₂), including the range of r₁₂ for which overlap occurs, must be incorporated into models of protein-protein interactions, and they determine the predicted values of B₂₂.

Previous efforts to model W and hence B₂₂ have been based on highly idealized descriptions of proteins. The protein molecules are almost always treated as spheres, although Vilker et al. (1981) also modeled bovine serum albumin (BSA) as an ellipsoid instead. For a sphere, the excluded volume is simply four times the sphere volume, but the appropriate sphere radius is uncertain. In addition, a steric barrier in the form of a "hydration layer" is sometimes assumed to surround the molecules and effectively represents an additional excluded volume; the thickness of the hydration layer may then serve as an adjustable parameter. The remaining contributions to B₂₂ are energetic, and they are usually modeled using standard colloidal methods (e.g., Hunter, 1987) to account for van der Waals (dispersion) attraction within the Lifshitz-Hamaker framework, and for electrostatic repulsion (Gallagher and Woodward, 1989; Muschol and Rosenberger, 1995). More elaborate treatments (e.g., Vilker et al., 1981; Haynes et al., 1992; Coen et al., 1995; Kuehner et al., 1997) may include such effects as permanent and induced dipole interactions, charge fluctuation terms, and osmotic attraction (depletion flocculation).

The differences among the various prior studies lie mainly in the assumptions used to calculate W; the idealized geometry is retained in all cases. Results are fitted to experimental data by balancing repulsive (mainly electrostatic) against attractive (mainly van der Waals) interactions, with the magnitude and/or the extent of one or more interactions treated as an adjustable parameter. Although these models have been partly successful in capturing some of the principal experimental trends and in obtaining reasonable numerical values with adjustment of parameters, discrepancies remain. For example, the Hamaker constant, which characterizes the magnitude of van der Waals interactions, can be used as an adjustable parameter, but fits to experimental data typically yield unrealistically high values. Among the reasons that have been identified for the discrepancies are the approximations inherent in describing the protein shape, in calculating the intermolecular interaction mechanisms, and in omitting effects that are poorly understood quantitatively, e.g., solvation interactions. A more subtle source of error is that W is usually found by adding independently determined contributions attributed to different mechanisms, each of which is found from the two-molecule pair potential by suitable orientational averaging. A more appropriate approach would be to add all of the contributions, accounting for their orientation dependence, before averaging.

The goal of the present work is to relax some of the assumptions previously made in computing B₂₂. Specifically, we account for details of shape and charge distribution by using crystallographic data. We have previously shown that a detailed shape representation gives rise to an excluded volume that is higher by ~70% than that for spheres of equivalent volume (Neal and Lenhoff, 1995), and that there is a pronounced orientational dependence for van der Waals interactions within the Hamaker framework (Roth et al., 1996). Here we extend those results by incorporating orientation-dependent interactions in calculating B₂₂. An indication that the orientation dependence may be important comes from an intriguing correlation between B₂₂ and protein crystallization, the latter of which is governed by close, specific contacts. Solution conditions conducive to crystallization have been shown to correspond to slightly negative values of B₂₂ (George and Wilson, 1994; George et al., 1997). Such conditions denote weak attraction, whereas stronger attraction (more negative B₂₂) was found to correlate with amorphous precipitation. Rosenbaum et al. (Rosenbaum et al., 1996; Rosenbaum and Zukoski, 1996) extended the correlation by showing that the measured virial coefficients could be used to predict solubilities within the framework of the sticky hard sphere model (Baxter, 1968), suggesting that short-range attraction is dominant in governing phase behavior. However, in view of the neglect of orientational dependence, such averaged interpretations may provide an incomplete picture of the underlying physics. This is indeed what emerges from our work, even though our models still rely on some obvious approximations because of both computational limitations and uncertainties regarding physical descriptions.

Two recent sets of virial coefficient measurements obtained in our laboratory by light scattering and small-angle neutron scattering (Velev et al., manuscript submitted for publication) provide a specific focus for our analysis. For hen egg lysozyme, B₂₂ decreases monotonically with pH (up to close to the pI) at each of a number of values of ionic strength between 0.005 M and 0.5 M, although the slopes of the curves become rather small at the higher ionic strengths. Within experimental error, the respective curves do not intersect. These results are qualitatively consistent with even the simplest models discussed above: both the decrease in charge with increasing pH and the increased screening with increasing ionic strength reduce repulsive electrostatic interactionshence the decrease in B₂₂.

For bovine chymotrypsinogen, there is again a monotonic decrease in B₂₂ between about pH 3 and 7, consistent with decreasing net charge, but in this case the curves at different ionic strengths all cross at around pH 5. Thus between this region and the pI (> 8), B₂₂ increases with increasing ionic strength, indicating that attractive electrostatics are dominant despite the finite positive charge. Similar behavior has been reported for bovine chymotrypsin (Coen et al., 1995), with model results attributing the trend to dipole interactions. However, the calculations described below, which are applied specifically to chymotrypsinogen, link the effect more closely to the geometric complementarity that is essential for specific molecular recognition, such as that manifested in protein crystals.

	THEORY AND METHODS

Top Abstract Introduction Methods Results & Discussion References

Calculation of B₂₂

(3)

(4)

Electrostatics

Electrostatic interactions can be quite realistically described and are essential for capturing the pH and ionic strength effects observed experimentally. The model used is a continuum one (Kirkwood, 1934; Warwicker and Watson, 1982; Sharp and Honig, 1990) in which each protein molecule is treated as a dielectric cavity (of dielectric constant 4) containing a set of point charges representing the ionized residues. Thus the charge distribution varies with pH. The molecules are immersed in an unbounded aqueous electrolyte solution. The interaction free energy, representing the electrostatic contribution to the PMF W, can be calculated once the potential field  throughout the system is known (Sharp and Honig, 1990).  is found by solving Poisson's equation within each protein molecule and the Poisson-Boltzmann equation (Hunter, 1987) in the electrolyte solution. Boundary conditions specifying continuity of the potential field and of the normal gradient of the electric displacement are satisfied at the dielectric boundaries.

The calculations performed here for chymotrypsinogen were based on chain A of the crystal structure in Protein Data Bank file 2CGA (Fig. 1). Charges were assigned to the ionizable residues according the Henderson-Hasselbalch equation by using intrinsic pK_a values (Zubay, 1984; Stryer, 1988). The resulting charges, 19.2e at pH 3 and 4.2e at pH 7, are within ~1e of the measured titration charge (Marini and Martin, 1971). The potential field was calculated by a boundary element method (BEM) (Yoon and Lenhoff, 1990), with the field in the presence of two molecules found iteratively (Zhou, 1993). For this purpose the molecular surface was discretized by using a cubic mesh (Zauhar, 1995), and the dependent variables (potential and normal flux) were approximated as varying linearly over each element (superparametric implementation (Aliabadi and Hall, 1988; Chan et al., 1990, 1992)).

View larger version (134K): [in this window] [in a new window]   FIGURE 1   Three-dimensional structure of bovine chymotrypsinogen A used as basis for calculations of second virial coefficient.

For calculations incorporating the full geometry of chymotrypsinogen, 2553 primary nodes, corresponding to 5102 cubic triangular elements, were used. The electrostatic computations for a system of this size would make the computational time prohibitive for the configurational exploration needed to evaluate the second virial coefficient via Eq. 4. Thus a more efficient approximate description was developed in which the spatial distribution of charges in each protein molecule was retained, but within a spherical dielectric cavity, the radius of which was based on the molecular volume of the protein. In this case, the charges were placed at the same angular positions in the sphere as in the full protein, with each charge at the same distance from the sphere surface as its distance from the nearest BEM mesh point in the full protein case. Because the tessellation used for the spheres comprised 362 primary nodes and 720 elements, the calculations using the sphere approximation were much more economical than those based on the full geometry. The approximation was thus used for the configurational integration in Eq. 4, but an extensive set of calculations using the full geometry was performed to aid in evaluating the accuracy of the approximation.

Short-range interactions

The second class of interactions includes a variety of short-range effects, such as van der Waals (dispersion), solvation, and hydrogen bonding interactions. These are collectively less well understood quantitatively than electrostatics, and because of the short range over which they act, they are very sensitive to the local geometry. These forces are significant in high-affinity protein-protein interactions such as antigen-antibody binding (Davies and Cohen, 1996), as well as in protein folding, because of the dense packing of the protein interior (Creighton, 1993).

Van der Waals interactions are the best understood of the short-range interactions. At the atomic level the interaction potential between atoms i and j (equivalent to the free energy for dispersion interactions) is usually described by a relation of Lennard-Jones form,

(5)

that captures the attractive nature of van der Waals interactions and the very short-range Born repulsion due to overlap of electron clouds. Here is the minimum interaction energy and _ij is the collision diameter. The drawback of relations of this kind for describing protein-protein interactions is that water molecules must be included explicitly, greatly complicating the computational task. The Lifshitz-Hamaker approach (Hunter, 1987) widely used in colloid science overcomes this difficulty via a continuum formulation in which the (usually attractive) interaction free energy of two bodies, in this case protein molecules, is described by

(6)

The integral is a function purely of the particle geometries, whereas the material properties are represented by A₁₂, the Hamaker constant, which is calculated from the optical properties of the interacting bodies and the intervening medium (water in this case). This approach also has a shortcoming, namely the breakdown of the continuum approximation at very short range. More generally, although Eqs. 5 and 6 are intended to describe the same physical behavior, the results of the two treatments are not easily reconciled.

Other short-range interactions present additional problems. Hydrogen bonding and solvation effects depend strongly on water and protein structures, and although atomistic simulations using Eq. 5 and pairwise Coulombic electrostatic models are quite successful at capturing local effects (Garde et al., 1997), straightforward means of extracting protein interaction energies are not yet available.

In the face of these difficulties, previous efforts to model the second virial coefficient have generally used Eq. 6, with the protein molecules treated as spheres. More detailed calculations (Roth et al., 1996) in which the idealization of shape is relaxed show that the geometric factor in Eq. 6 is almost always overestimated by the sphere model relative to the value for a more realistic geometric representation. Even here, though, questions remain because of the breakdown of Eq. 6 at short range and the neglect of other interactions. Because of these difficulties, interaction energies for high-affinity contacts are expediently calculated by empirical approaches based on some measure of the complementarity of the apposing surfaces. The buried surface area is frequently used (Novotny et al., 1989; Horton and Lewis, 1992), based on the correlation observed between the accessible surface area of functional groups (usually nonpolar) and the free energy of transfer between water and an organic solvent (Creighton, 1993). An alternative measure of complementarity is the volume of the intervening "gap" or a "gap index," calculated as the ratio of the gap volume to the buried area (Jones and Thornton, 1996).

Because we need to calculate the interaction energy for two molecules not only in a bound state, but also in arbitrary configurations, we use a different semiempirical approach for describing the short-range interactions. First, where the intermolecular gap is large enough to justify a continuum description, we use the Lifshitz-Hamaker approach as applied to realistic representations of the protein shape (Roth et al., 1996), with a protein-water-protein Hamaker constant calculated to be 3.1kT (Roth et al., 1996). At shorter range, however, we use Eq. 5 to capture short-range attraction and repulsion. Even in this situation, though, substantial parts of the two protein molecules are widely separated. Thus our method is a hybrid one in which all atom pairs separated by more than an arbitrary center-to-center distance (here chosen to be 6 Å, the approximate distance that allows a water molecule to fit into the gap) are assumed to interact by Eq. 6, and all other pairs by Eq. 5. Although the transition is taken to be discontinuous, the protein-protein interaction energy curve is fairly smooth as a function of distance because of the large number of atom pairs involved.

This formulation is clearly far from rigorous, but it captures surface complementarity very well and provides a fairly smooth interaction energy profile for two molecules being translated toward each other. A specific advantage over the usual Hamaker approach is that there is no singularity at contact; instead, there is a smooth minimum in the energy trajectory, and for several crystal contacts and high-affinity contacts that we have examined, this minimum lies within a few tenths of an angstrom of the crystallographically observed contact (Neal, 1997).

We have used the OPLS parameter set (Jorgensen and Tirado-Rives, 1988) as the basis for our short-range calculations with Eq. 5, but based on measured binding energies of high-affinity pairs (e.g., Bhat et al., 1994), the resulting interaction energies are too large in magnitude. This discrepancy is not surprising given the important effects, notably those involving solvation discussed above, that are neglected in our short-range calculations. We find empirically that a reduction in the contribution of Eq. 5 by a factor of 2 reduces the short-range energy to a more realistic value; interestingly, this rescaling brings the Lifshitz-Hamaker and OPLS energies calculated in vacuum (Roth et al., 1996) into good agreement.

	RESULTS AND DISCUSSION

Top Abstract Introduction Methods Results & Discussion References

Because of the emphasis in our calculations on seeking to explain experimental results apparently associated with attractive electrostatic interactions (Coen et al., 1995; Velev et al., manuscript submitted for publication), we first present results of the electrostatics computations and then examine how they contribute to B'₂₂ via application of Eq. 4. Electrostatics calculations using the complete description of protein shape and charge distribution were performed for chymotrypsinogen at 0.1 M ionic strength for pH 3 and 7. For each set of solution conditions, five separation distances at each of 36 angular orientations were sampled; the separation distance was based on the minimum distance between the van der Waals surfaces of the nearest atom pairs. The ranges of the electrostatic interaction energies obtained are shown in Fig. 2. At pH 3 all interactions are repulsive, as would be expected from the almost complete protonation of the acidic residues. The interaction energies cover a fairly wide range, however, which is attributable to differences in the proximity of positive charges on each pair of molecules in the different configurations. At pH 7 there is a qualitative difference: in at least one configuration the interaction is attractive, despite the positive net charge. An interesting feature is that the magnitudes of the repulsive interaction energies are similar at the two pH values, despite the large differences in net charge. This and the appearance of attractive interactions reflect the importance of the more closely spaced pairs of individual charges as opposed to the global net charge.

View larger version (29K): [in this window] [in a new window]   FIGURE 2   Range of calculated electrostatic interaction energies for chymotrypsinogen as a function of separation distance at pH 3 and pH 7, 0.1 M. Full protein structural details were included. Inset shows distribution of energies at 1-Å separation for the 36 configurations studied at each pH value. Filled bars, pH 3; unfilled bars, pH 7.

The corresponding calculations for the spherical geometry with distributed charge were initially performed for a sphere of volume equal to that of the protein molecule (19.4 Å radius). The approximate model was found to overpredict the electrostatic interaction energy in all but one orientation, by as much as a factor of 2. Imoto (1983) suggested that, in such an approximate model, the spacing between the charges should be kept as close to that in the true protein structure as possible. To achieve this without placing charges outside the low dielectric cavity as Imoto did, we used a larger sphere radius. Calculations of the excluded volume of protein molecules (Neal and Lenhoff, 1995) indicate that from the point of view of protein-protein interactions, the cavity occupied by the protein has a radius ~20% greater than that of the sphere of equivalent volume. Thus the sphere radius was increased to 23.1 Å while keeping the charge depth below the surface and the angular distribution fixed. The results for pH 7 (Fig. 3) show that the sphere model still overestimates the electrostatic energy, but by a smaller factor; comparable results were obtained for pH 3. Although it may be possible to improve agreement further by manipulating the sphere radius, the algorithm used to generate Fig. 3 was used for the full computations because it represents the orientational dependence quite well, and the actual discrepancies are fairly small given the other approximations made, especially in the short-range interaction calculations.

View larger version (23K): [in this window] [in a new window]   FIGURE 3   Relation of electrostatic interaction energies calculated for full chymotrypsinogen representation (Fig. 1) and sphere approximation (radius 23.1 Å). pH 7, 0.1 M.

The electrostatic attraction apparent at pH 7 in Figs. 2 and 3 would not be expected to dominate B'₂₂ behavior in terms of a simple average of the energy, but its importance is boosted by the Boltzmann weighting in Eqs. 3 and 4. However, this weighting applies to the total PMF W, and not to the individual contributions, so it is necessary to examine the short-range interactions as well. The critical attribute of these interactions is that they are attractive, except for the Born repulsion (see Eq. 5) at very short range, so that the integrands in Eqs. 3 and 4 should be positive if only the short-range interactions are considered. The values of I_in are then positive for all orientations, corresponding to negative contributions to B'₂₂. The squares in Fig. 4 show the calculated values of I_in, based on short-range interactions only, as a function of the well depth of the short-range interaction energy profile with r₁₂; 1792 configurations, sampled by a Monte Carlo integration routine, are shown. All of the points lie roughly on a straight line on this semilog plot, indicating that the value of I_in is controlled mainly by the short-range energy well depth rather than reflecting the full short-range energy profile; this is a consequence of the exponential Boltzmann weighting.

View larger version (36K): [in this window] [in a new window]   FIGURE 4   Relation of calculated inner integral to well depth of short-range interaction energy profile in 1792 angular configurations of chymotrypsinogen pairs. , Short-range interactions only; +, short-range plus electrostatic (pH 7, 0.1 M) interactions.

When electrostatic interactions are included in the B'₂₂ calculations as well, the I_in points in Fig. 4 are perturbed from the values previously examined. For a given configuration, attractive electrostatics increase I_in and repulsive electrostatics decrease it, giving rise to the points shown as + in Fig. 4 for electrostatics calculated at pH 7 and an ionic strength of 0.1 M. The extent to which each I_in changes depends on the relative magnitudes of the short-range and electrostatic energies. The values in the vicinity of the short-range well are shown in Fig. 5 for the different orientations sampled; there is clearly no correlation between the two energies. Electrostatic interactions are predominantly repulsive, and for some orientations they overcome the short-range attraction, giving rise to negative values of I_in; these are not shown in Fig. 4 because of the logarithmic scale used, but they are relatively small because of the Boltzmann weighting. The attractive electrostatic interactions, on the other hand, are rather weak (Fig. 5), so very large values of I_in are due mainly to strong short-range interactions. However, because of the Boltzmann weighting, attractive electrostatics can boost I_in substantially above the value due to short-range interactions alone (Fig. 4); as discussed below, these configurations may be central to determining trends in B₂₂.

View larger version (32K): [in this window] [in a new window]   FIGURE 5   Distribution of short-range interaction energy well depths and electrostatic interaction energies (pH 7, 0.1 M) in corresponding configurations.

An additional feature of the I_in values in Fig. 4 is the very wide range covered, which is a critical factor in determining B'₂₂. Equation 4 shows that B'₂₂ comprises the (positive) excluded volume contribution and an unweighted orientational average of I_in, and in view of the many orders of magnitude spanned by I_in, it is the orientations with the largest values of I_in that most strongly affect B'₂₂. This can be seen more clearly from the histogram in Fig. 6, showing the distribution of I_in values for the case shown in Fig. 4. Orientations for which I_in is smaller in magnitude than ~10⁵ ų contribute negligibly to the energetic term in Eq. 4; this includes those orientations for which repulsive electrostatics lead to negative values of I_in. Conversely, the orientations that are most important in determining B'₂₂ are those for which I_in is large, and as is apparent from Fig. 4, the main effect causing this is strong short-range interactions. Such interactions result from a high degree of complementarity between the apposing surfaces, characteristic of a molecular recognition phenomenon. Qualitatively similar results would be obtained even if the short-range interactions were calculated by some other approach, e.g., one based on buried surface area. What is essential here is to account in detail for the geometry of the interacting molecules, which is overlooked in models treating the protein molecules as spheres.

View larger version (37K): [in this window] [in a new window]   FIGURE 6   Distributions of inner integral values for 1792 angular configurations of chymotrypsinogen pairs, with different contributions included as shown in legend. Inset shows enlarged view of high-Iin region.

The observed effect on B'₂₂ of electrostatic interactions must then also be related to the high-complementarity configurations: it is the nature of the electrostatic interactions in these orientations that plays the main part in determining experimental electrostatic effects. Attractive electrostatics will increase the I_in from those obtained for short-range interactions only, leading to a decrease in the value of B'₂₂, and vice versa for repulsive electrostatics. Because of the nonlinearity of the Boltzmann factor, however, the effects of attraction and repulsion are asymmetrical, and the net result of the orientational average is not easily predicted. As is apparent from Figs. 4 and 5, however, in this system there is a disproportionate number of high-complementarity configurations for which the electrostatics are attractive at pH 7 and 0.1 M ionic strength. Such attractions are screened at higher ionic strengths, and although repulsive interactions are screened as well, the asymmetry of the Boltzmann dependence, together with the distribution and coupling of electrostatic and short-range energies (see, e.g., Fig. 5), can result in the experimentally observed increase in B₂₂ with ionic strength for chymotrypsin and chymotrypsinogen (Coen et al., 1995; Velev et al., manuscript submitted for publication). The histograms of I_in values in Fig. 6, all for the same set of orientations, show that there are more configurations with I_in > 10⁸ ų at 0.1 than at 0.3 M.

Beyond the above observations, B₂₂ cannot be estimated accurately, because of the uncertainties in the calculated interaction energies and the computational challenge of fully sampling the configurational space. However, the best results obtained using the Monte Carlo integration routine are in reasonable quantitative agreement with experimental results (Table 1), especially in view of the fact that no adjustable parameters were used in the calculations. The counterintuitive trend with electrolyte concentration at pH 7 is also correctly captured. Given the large error estimates on the calculated B₂₂, however, the agreement is significant not so much in the correct prediction of values and trends, but rather in the plausibility of the formulation used. The central role of short-range interactions that emerges from our work was postulated previously (Rosenbaum et al., 1996), but the more specific importance of the high-complementarity configurations represents a new paradigm for interpreting B₂₂ data. Another implication is that in view of the physical uncertainties, computational complexity, and dependence on detailed structural information required to calculate B₂₂ within this framework, detailed prediction of solution thermodynamic properties is extremely difficult, and one must rely on experimental probes of protein interactions.

The new paradigm also suggests a natural link that may explain the correlation between B₂₂ measurements and protein crystallization (George and Wilson, 1994; George et al., 1997): both are governed by molecular recognition in a relatively small number of configurations. To explore this further, we have compared the energetics for the high-complementarity configurations identified in the Monte Carlo integration (Figs. 4 and 5) with the crystal contacts in the crystallographic coordinate files (Wang et al., 1985). The most attractive crystal contacts yield short-range interaction energy well depths in the range 13-24kT. The corresponding values in a new crystal form of chymotrypsinogen that we have recently identified (Pjura et al., manuscript in preparation) are also on the order of 20kT. Thus although relatively few high-complementarity configurations have been been found in the configurational exploration, these do indeed have affinities that are comparable to those of crystal contacts. That there are more such configurations than actually observed in a given crystal structure is consistent with previous analyses (Janin and Rodier, 1995) and with the ubiquity of crystal contact configurations observed in the different crystal forms of some proteins (Crosio et al., 1992). Which ones are manifested in a given structure may be determined by such effects as kinetics (not considered here); electrostatic repulsion may play a role.

In summary, our results indicate that B₂₂ is determined largely by the contributions of relatively few, highly attractive configurations. The effects of varying electrostatic parameters such as pH and ionic strength are also manifested mainly because of the influence of these configurations. This central role of molecular recognition suggests a plausible link between solution thermodynamic properties and the formation of crystalline phases or other self-organized structures. Our conclusions also suggest that local structural perturbations are an underused route to manipulating protein interactions; these may be effected by, for example, specific adsorption of ions or surfactants. The results also indicate the vast complexity of molecular interactions that may occur in solutions of protein mixtures, not least the intracellular environment.