CRYSTAL GROWTH Journal of Crystal Growth 143 (1994) 71—78
ELSEVIER
Molecular interactions in crystallizing lysozyme solutions studied by photon correlation spectroscopy Wolfram Eberstein, Yannis Georgalis ~ Wolfram Saenger Institutfür Kristallographie, Freie Universitdt Berlin, Takustrasse 6, D-14195 Berlin, Germany Received 7 May 1994; manuscript received in final form 6 June 1994
Abstract The concentration dependence of the diffusion coefficient of particles suspended in solution depends primarily on the occupied volume fraction and on repulsive and attractive forces. This dependency is expressed by the interaction parameter that can be assessed experimentally by light scattering measurements. In the present work we have determined the interaction parameter for the diffusion coefficient of hen egg white lysozyme under nonstationary crystallization conditions. The hydrated lysozyme monomer has a free-particle diffusion coefficient D 2 s~that corresponds to a hydrodynamic radius a 2.09 nm. It also bears a net charge of +6.40 and102has±4a ~tm Hamaker constant, AH 7.7 kBT, for the crystallization conditions employed. At an intermediate degree of screening, Van der Waals forces become dominant, and the behavior of A can be correlated with the onset of fractal cluster formation. The latter is employed as a useful index for predicting the fate of the lysozyme solution under examination. =
=
=
1. Introduction In supersaturated protein solutions, monomers assemble to form aggregates that form either crystals or amorphous precipitate. Light scattering studies on the clusters formed during protein crystallization is currently an interesting topic under investigation in many laboratories. We have used hen egg white lysozyme as a model protein for investigating the conditions that determine crystal formation. We find that the interaction parameter A, which defines the concentration dependence of the diffusion coefficient (for a review, see Ref. [11),can be employed as an alternative way of diagnosing crystallization in the following sense: At moderate volume frac_______
*
Corresponding author,
tions, unless the repulsive electrostatic potential is extinguished and the attractive Van der Waals forces become effective, the slowness of the rate limiting step (dimerization) precludes the onset of crystallization. This condition may be necessary but it is not sufficient for crystal growth and guarantees only the onset of aggregation. Whether such aggregates will lead to crystal, precipitate or gel has been recently investigated by our group [2—61.
2. Theory 2.1. Photon correlation spectroscopy (PCS) The theoretical aspects of the technique as well as the data evaluation algorithms employed
0022-0248/94/$07.OO © 1994 Elsevier Science B.V. All rights reserved SSDI 0022-0248(94)00359-T
72
W. Eberstein et al. /Journal of Crystal Growth 143 (1994) 71—78
in the present work have been reviewed recently by Brown [71and need not be discussed again. PCS provides a direct means for determining the diffusion constant of an infinetely dilute suspended particles through the rate constant F of the decay of the scattered field photon correlation function 0,
kBT/6~a.
=
=
(2)
Here kBT denotes the thermal energy and ~ the viscosity of the solvent.
4f[g(r) a o
2.2. First order interaction parameter For concentrated systems these simplified assumptions are not valid, and the various types of interactions between particles have to be taken into account. The problem of interactions, as determined by PCS, has been under investigation for over twenty years, but theories disagree in respect of both the magnitude and the sign of A, as shown in a compilation given in Ref. [8]. Through the years, the most viable theories are the ones developed by Batchelor [9,10] and Felderhof [11,121. For concentrated macromolecular solutions in the hydrodynamic regime, Eq. (2) has to be modified to account for finite concentration effects, correct to first order in concentration D°~ D0[1 =
D0[1
+
A4S
+ o(~2) +
. .
+ (A~+ A0 + AA + 2) + + O(4
]
—
hr2 dr,
(4)
3
AA=
As AD
=
=
—
3 ~19 a ~f [~(7) a
1]r dr,
6 -
(5)
5 a ~(7) }g(r)r2
dr,
(6)
~a4f ~dr, or 1.
(7) (8)
For hard, non-interacting spheres g(r) 0 if r < 2a or g(r) 1 if r ~ 2a and Ày 8, A0 AA 1.73 and A~ 0.29. These contributions =
=
=
=
—
A
(1)
where q denotes the scattering vector and D0 the translational diffusion coefficient that can be directly associated with the hydrodynamic radius a of the free particle via the relation D0
A~
0= —~-f[g(r) a o
2D F=q
be individually expressed as functions of the density distribution function g(r) [13,141as follows
—
=
=
=
and AD sum up to a constant value of 1.56 that accounts only for the excluded volume effect. An equivalent formulation has been given by Corti and Degiorgio [15] using the results obtained by Batchelor [9]. The density distribution function g(r) is related to the pair potential of mean force V(r) [16,17] via g(r)
=
(9)
~
therefore, if V(r) can be realistically estimated, the potentials involved can be compared with computer simulations of three body correlations. Theoretical developments made along these lines of research have been reviewed recently in Refs. [18,19]. 2.3. Energy of interaction
.1 A~+ AD)~ (3)
For calculating the interaction parameter A, the effective potential V(this r) should be known. As ainteraction first approximation potential can be considered, according to the DLVO the-
In Eq. (3), A~and A 0 denote the virial and the Oseen parts, assuming stick boundary conditions for the hydrodynamic interactions, AA and A~ denote the self-interaction and short-range parts, and AD the dipole part. These contributions can
ory [201,as the sum of two contributions involving the attractive Van der Waals energy VA and the repulsive Coulombic energy VR, expressed as a function of the dimensionless quantity x (r 2a)/2a. =
—
W. Eberstein et al. /Journal of Crystal Growth 143 (1994) 71—78
The attractive forces acting between two molecules in solution are caused by the Van der Waals induced dipole interaction: AHI
VA(x)=
1
—~[X2+2X
+ 2 in
+
(x
4. Materials and methods (10)
+ 1)2
where AH denotes the Hamaker constant [21]. Eq. (10) is divergent for x 0 and the numerical implementation of this expression requires a lower cut-off value that may considerably modif3T the magnitude of the attractive energy VA(x). We have increased the lower integration bound by a value that could approximately qualify for the thickness of one Stern layer. At low particle and ion concentrations, the potential of mean force can be approximated by the well known mean field expression =
VR( x)
=
where
K
K
=
(
~
2
ln( 1
+
e
—
(11)
2Kax),
is the inverse Debye—HUckel length [21],
e2NA LZ2n J
1/2 JO)
increase of the particle radius due to the bound counterions. Z~,denotes the charge of the protein, and a5 the radius of the screening countenon that is taken as 0.18 nm for Cl.
1 2 (x+1)
\1
(x~+2x
73
(12)
,
4.1. Materials
All chemicals used were of analytical grade. Water was obtained from a Millipore-Q device. Hen egg white lysozyme, three times crystallized, was purchased from Sigma Chemical Co. (Deisenhofen, Germany). Protein preparations were treated as previously described [2—6]and a standard 0.1M Na-acetate buffer, adjusted to pH 4.2, was employed throughout. Before aggregation experiments were initiated, the preparations were tested for monodispersity with PCS in the absence of salt and were found to be monodisperse at the highest protein molanity employed. Equal volumes of salt and protein solutions, at twice the desired molarity, were rapidly mixed and filtered through Minisart (Sartorius GmbH, Germany) sterile filters, 0.22 ~m pore size, into standard cylindrical light scattering cells. 4.2. PCS data acquisition and evaluation —
and e denotes the electronic charge, NA Avogadro’s number, the permittivity of the solvent,
PCS measurements were conducted with an
Z
1 the ion valence and n10 is the number concentration of the Jth ion. We have accounted for the small drop of the permittivity with added electrolyte by using a simple empirical correction e~(1 0.75c5) [22], where ~ is the dielectric constant at zero salt content and c~the salt molarity. For a given protein charge, an increase of the salt molarity increases the value of K. Because protein charge depends also on pH, salt concentration is not the only variable that influences the magnitude of K. In Eq. (11), ~ denotes the surface potential, given by =
—
ZP e e
(13 )
The diffusion coefficient of the monomers was calculated byCONTIN Laplace-inverting the spectra with the program [23—25]implemented in a
The last factor, in brackets, accounts for the
CONVEX-C220 supercomputer. CONTIN re-
ip
0
e
F a
I
e’~° +as) + K(a + a5)
____
1
1
ALV/SP-86 spectrogoniometer (ALV, Langen, Germany) and the ALV-FAST/5000 digital autocorrelator boards at a scattering angle of 20°and a constant temperature of 20 ±0.1°Cthroughout. A tunable Ar~laser (488 nm) was employed as light source. For each experiment thirty measurements of 15 s duration each were performed. Experiments were initiated within less than 1 mm after mixing the solutions. This was necessary since under the employed conditions large fractal aggregates, that dominate the scattering process, appear within a few minutes after mixing protein and salt.
74
W. Eberstein et al. /Journal of Crystal Growth 143 (1994) 71—78 [mMI
C p 00
1.0
~
~.o
-
3.0-
r
H 5,0
1 .20 I
- -~
1. If we set the correlation time F~ equal to X T, 2 2.7 for our conditions we obtain ‘r 1/D0q i0~ s. Following Ref. [1], the magnitude of TB can be calculated from the mass of the lysozyme monomer Mm 14600 and the free particle friction coefficient ~ 6ri,a as TB Mm/to 6.33 x iO~ s. The upper limit for T 1 can be taken as the time required for a lysozyme molecule to diffuse the average distance between molecules, say in the most dilute protein solution used in our =
~“4.0
.,---‘~
~-~‘~H-
-
=
=
=
1.05
0.00
0.02
0 (14
0,08
6 0.08
8
=
=
x 10—6 cm. experiments, i.e., 3=1.5 c~ 0.5mM. ThisTherefore distance r~ is can calculated as r2/6D thenber=Kc~Y~ 0 3.97 x iO~s, and it follows that T>> T1 >> TB. More elaborate treatments [27—301have been developed for the computation of the total interaction energy V(x) and for further assessment of =
0.10
=
Fig. 1. Normalized diffusion coefficient of lysozyme monomers as a function of the protein volume fraction 4 (lower x-axis) or protein molarity c~ (upper x-axis). Each point is the average of at least twenty measurements. Error bars denote three standard deviations. At later times, i.e. 20 mm, the scattering amplitude of the monomers is very low and renders reliable inversions difficult. The slopes of the curves yield the interaction parameter through Eq. (3); NaC1 molarities are: 0.05M (., 1); 0.1M (v, 2); 0.2M (v, 3); 0.3M (0, 4); O.5M (E, 5); 0.8M (~,6); 1.1M (~,7) and 1.4M (~,8).
solves well the monomeric species during the first few minutes of the experiment even when high salt and high protein molanities are employed. Eight different NaC1 molanities in the range between 0.05M and 1.4M have been studied. Experiments conducted at constant NaCl molarities involved six different lysozyme molanities in the range between 0.5mM and 3.0mM; results are displayed in Fig. 1. The interaction parameter A was obtained by a simple linear regression of D~ versus 4 for each NaC1 molarity. Higher order fitsinfinite were not justified. respective at dilution haveThe been employedintercepts for normalizing the data to the dimensionless ratio
D~/DO.
the static structure factor or interaction parameter. However, the inherent uncertainty of our non-stationary conditions, i.e. determination of the cooperative diffusion coefficient in the presence of large aggregates, at present does not justify a higher level of complexity. We reckon that a simple model like the one adopted above will be sufficient for capturing the elementary physics of the problem. 5.1. Monomer diffusion coefficient D0 and interaction parameter A The inverse Debye—Hückel screening length the interaction parameter A, and the monomer diffusion coefficient, extrapolated at infinite dilution, D0, are given in Table 1. The value of the free-particle diffusion to coeffi2 s’ extrapolated zero cient,molanity D0 102and ±4 zero ~m volume fraction, agrees salt with the values reported by other investigators [31—33], 103—106 ~~m2 s~ for the lysozyme monomer. The hydrodynamic radius resulting from Eq. (2) is a 2.09 nm. We can now compare this value of the hydrodynamic radius with the minimum radius expected for the equivalent unhydrated sphere K,
=
=
5. Results The experiments were conducted in the longtime regime [26] that is characterized by time scales longer than the Brownian motion time scales of the particles T>> ‘TB’ and longer than those of the interactions between particles, T>>
am 3. (14) 10 (3Mm3/4~rNA)~ Using a partial specific volume of i~ 0.703 cm3 g~ [33] and a mass of 14600 g for the lysozyme =
=
W Eberstein et al. /Journal of Crystal Growth 143 (1994) 71—78 Table 1 Inverse Debye—Huckel screening length ic~interaction parameter A, and infinite dilution extrapolated diffusion coefficient D 0, for eight NaCI molarities; the interaction parameter has been extracted from linear regression of measurements sets involving at least six different lysozyme concentrations in the range between 0.5mM and 3mM to allow an accurate determination to be made of the collective diffusion coefficient Dc as a function of volume fraction ~ CNaCI
K
(M) 0.05 0.10 0.20 0.30 0.50 0.80 1.10 1.40
(i0~m’) 0.73 1.04 1.46 1.79 2.31 2.93 3.43 3.87
A
D0 2/s) (~m 99.6±2.7
+2.55±0.81 +1.76±0.61 — 1.67±0.35 —2.89±0.57 —3.63±0.34 —3.69±0.40 3.71±0.71 —3.82±0.56
97.1±2.2 104.0±1.4 98.3±3.9 102.0±1.4 94.3±1.7 97.4±1.9 93.1±4.4
monomer, we obtain amjn 1.59 nm. The ratio between the hydrodynamic radius- and the radius of the equivalent sphere is 1.3. The crystal structure shows lysozyme to be an ellipsoid of revolution with semiaxes a 2.25 nm, f3 y 1.50 nm; we can calculate an equivalent hydrodynamic radius to be 1.79 nm. If we take the thickness of one water layer to be 0.28 nm, we obtain a value of 2.07 nm, which agrees closely with the determined value of 2.09 nm. =
=
=
=
75
identified. Therefore a net charge of + 7 is expected. The curves displayed in Fig. 2 intersect at a common locus with a protein charge Z~ + 6.4 and a Hamaker constant of A11 7.7 k BT. Slightly different values + 6.2 and AH 8.4 kBT were obtained when using the Batchelor— Felderhof results in the form given by Corti and Degiorgio [15]. The value determined for the charge agrees with the net charge predicted from the crystal structure if one considers the approximations and assumptions involved. The magnitude of the Hamaker constant is also close to the values reported for aqueous polystyrene solution [36,37], 0.26 kBT to 5.7 kBT, those of ionic micelles 11.3 kBT [15], and those determined for hydrocarbons =
=
=
=
and their derivatives [38,39]. For the determined value of Z~,we can calculate the surface potential, ~ to vary between 15 and 27 mV (Eq. (11)) which are physically reasonable values. The deviation observed for the data at ~a 1.45 results from the low value of ~a. This value is close to one where the repulsive potential is ill-defined and Eq. (11) is not immediately applicable. Further, it is also implicitly assumed, as a very crude approximation, that the electrostatic potential is centrosymmetric. Finally, at low salt concentrations, minor populations of dimers or =
5.2. Molecular parameters z and AH
____________
25
__________________ -—
-
/1
The interaction parameter A derived from Eqs. (4)—(7) is a non-linear function of the Hamaker constant A 11 and of the charge Z~of the protein. To obtain the correct values of these parameters, matrices of A values were computed for different ionic strengths and for wide ranges of A11 and Z~.Computations of the integrals were made with a Romberg quadrature [34] in a CONVEXC220 supercomputer. Pairs of the values (x =AH, y Z~)satisfying Eqs. (4)—(7) can be found by calculating the intersection of the x, y surface with z A. These pairs of values for different salt molarities are shown in Fig. 2. From the known crystal structure of lysozyme [35] (Protein Data Bank entry 2LYZ), 16 positive and 9 negative charges have been =
=
20
15 10
0 0
2
4
6
8
10
12
14
16
z I’ Fig. 2. Hamaker constant A11 versus protein charge Z~,.The plotted curves correspond to different the measured values of the interaction parameter A at the degrees of screening indicated in Table 1. The curves have a common locus with coordinates Z~,= +6.4 and AH = 7.7 kBT.
76
W. Eberstein et al. /Journal of Crystal Growth 143 (1994) 71—78 4
_______________________
-- ______________
__________
_
2
____________
~j69~p~-
0 —1
S.
—2
_
—4 —3 —5
1
2
H’ 3
4
5
6
7
8
io2 i02
9
t
~ a
Fig. 3. Interaction parameter as a function of salt molarity. The error bars indicate three standard deviations and the solid line was calculated using Eqs. (4)—(7). Charge Z~,= + 6.4 and Hamaker constant AH = 7.7 kBT.
other higher oligomers can cause sufficient errors in the scattered intensities. Essentially, this happens due to large excluded volume effects, because the screening length of the interactions increases with decreasing salt content. Therefore such deviations are not unexpected at lower screening lengths. The interaction parameter, A, is displayed in Fig. 3 as a function of the dimensionless product ~a. The solid line is calculated with the above determined values of protein charge and Hamaker constant and satisfies the observation. It drops
0.8
1
~
2
0.4 0.0
[~]
Fig. 5. Typical aggregation kinetics of 2.1mM lysozyme with 0.75M NaC1, pH 4.2, plotted in a double logarithmic fashion. 1”j following the The solid line is a power-law fit R(t) ~ t procedure described in Ref. [21.The fractal dimension determined after decoupling the rotational motions is df = 1.69 and the minimum (zero time extrapolated) cluster radius is R 0 = 12 nm.
from positive to negative values and reaches a constant value at about ia 4.6 that corresponds to a concentration of 0.5M NaC1. The calculation of molecular data allows the evaluation of both repulsive and attractive potentials as a function of screening, Fig. 4, within the limitations of the theoretical approximations used. At low NaCI molarity, 0.05M, the electrostatic repulsion poses an activation barrier to dimerization of only a fraction of the thermal energy 0.5 kJ/mol. Above 0.5M NaCI, screening is virtually complete. Below 0.5M NaCI, the examined solutions remain stable (only minor populations of smaller clusters are detectable) but crystal growth is not observed. In the salt range of 0.5M—1.OM crystals are formed within one or two days, and =
above 1 .5M the protein precipitates irreversibly after a few hours. Results from detailed kinetic experiments have been given elsewhere [2—6]and a typical kinetic experiment recorded during these measurements is displayed in Fig. 5.
—1.6 —2.0 0.00
0.25
0.50
0.75
1.00
1.25
6. Conclusions
x
Fig. 4. Calculated interaction potentials of lysozyme. NaC1 molarities are: 0.05M (1), 0.1M (2), 0.2M (3), 0.3M (4), 0.5M (5), 0.8M (6), 1.1M (7), and 1.4M (8).
The analysis of light scattering data on systems of interacting particles can only be carried out meaningfully if care is exercised in the computa-
W Eberstein et al. /Journal of Crystal Growth 143 (1994) 71—78
tion of interactions. For charged proteins, these interactions are dominated by Coulombic repulsion. As the ionic strength is increased, Coulombic interactions decrease and Van der Waals interactions prevail; particles collide with each other and build dimers or higher aggregates. This effect manifests itself as an increment of the scattered intensity and as a decrease in the collective diffusion coefficient. Under optimal crystallization conditions, energetic barriers should be as low as possible to allow dimerization and further nudeation to take place. In the present work we have adopted a simple treatment for interpreting the PCS results and obtaining estimates of the charge and of the Hamaker constant of lysozyme. This is a very important step toward the diagnostics of protein crystallization, since we can assess the instability threshold for a given precipitating agent. While Van der Waals interactions do not influence the diffusion coefficient of stable suspensions during nucleation, the system is far from equilibrium and it is expected that Van der Waals forces will be operative when Coulombic forces expire. In the limit of infinite dilution, at pH 4.2, lysozyme behaves as a hydrated monomer with diffusion coefficient of D0 102.3 ~tm2 s’, bears a net charge Z~, + 6.4 and a Hamaker constant A 11 7.7 kBT. Both values fit well the experimentally determined interaction parameter A when the latter is determined as a function of ~a. These molecular parameters will be useful for computer simulations concerning the aggregation and crystal formation processes. While for low molecular weight systems crystal growth can be theoretically described and to some extent predicted and manipulated, these interactions are far too complex for systems of biological macromolecules and only poorly understood. For the lysozyme—NaC1 systems the interaction parameter A can serve as a potential index of the initiation of aggregate formation since it signals the expiration of repulsive forces. The formation of fractal structures has often been questioned in batch crystallization experiments. In this work we present persuasive evidence that these structures are not artifacts and put the reasoning for their formation in a sound frame. =
=
77
The DLVO framework, used here to model the interactions, should be considered only as a rough, first order approximation. It neglects both non-Coulombic ion—ion repulsion and ion—ion correlations, and overestimates the importance of the Van der Waals interactions, since it is assumed that they are not modified by the solvent. A more elaborate treatment of the problem incorporating recent advances in the theory and simulation of electrolytes [40] will be presented in a subsequent communication (Soumpasis and Georgalis [41]). We have shown that the concentration dependence of the diffusion coefficient of monomeric species in aggregation experiments under optimal or suboptimal crystallization conditions may be used as a primary tool for obtaining information concerning the future of the solution. At least in experiments with lysozyme there is a strong correlation between the determined interaction parameter, A, and the ability of the solutions to produce well diffracting crystals. These studies have to be extended to other proteins that serve as candidates for studying protein crystallization. Acknowledgements
=
We wish to thank Drs. M.D. Soumpasis, D.N. Petsev, P.R. Wills and P. Zielenkiewicz for advice and criticism on the manuscript, and for making available unpublished works. We also thank the European Community and the European Space Agency (ESA/ESTEC) for financial support.
References [1] P.N. Pusey and RJ.A. Tough, in: Dynamic Light Scattering, Ed. R. Pecora (Plenum, New York, 1985) p. 85. [2] Y. Georgalis, A. Zouni, W. Eberstein and W. Saenger J. Crystal Growth 126 (1993) 245. [3] Y. Georgalis and W. Saenger, Adv. Colloid Interf. Sci. 46 (1993) 165. [4] W. Heinrichs, Y. Georgalis, H-i. Schbnert and W. Saenger, J. Crystal Growth 133 (1993) 196. [5] W. Eberstein, Y. Georgalis and W. Saenger, Eur. Biophys. J. 22 (1993) 359 [6] Y. Georgalis, J. Schüler, W. Eberstein and W. Saenger, in: Fractals in the Natural and Applied Sciences, Ed.
78
W Eberstein et al. /Journal of Crystal Growth 143 (1994) 71—78
M.M. Novak (North-Holland—Elsevier, Amsterdam, 1994) p. 139. [7] W. Brown, Ed., Dynamic Light Scattering, The Method and Some Applications (Oxford Science Publ., 1993). [8] T.G.M. van de Ven, Colloidal Hydrodynamics (Academic Press, New York, 1989). [9] G.K. Batchelor, J. Fluid Mech. 74 (1976) 1. [10] G.K. Batchelor, J. Fluid Mech. 131 (1983) 155. [11] B.U. Felderhof, Physica A 89 (1977) 373. [121 B.U. Felderhof, J. Phys. A (Math. Gen.) 11(1978) 929. [13] J.M. Ziman, Models of Disorder (Cambridge University Press, 1978). [14] R.J. Hunter, Foundations of Colloidal Hydrodynamics (Clarendon, Oxford, UK (1989). [15] M. Corti and V. Degiorgio, J. Phys. Chem. 85 (1981) 711. [16] J.G. Kirkwood, J. Chem. Phys. 2 (1934) 351. [17] J.G. Kirkwood, Chem. Rev. 19 (1936) 275. [18] G. Hummer and M.D. Soumpasis, Mob. Phys. 75(3) (1992) 633. [19] G. Hummer and M.D. Soumpasis, J. Chem. Phys. 98 (1993) 581. [20] E.J.W. Verwey and J. Th. G. Overbeek, Theory of the Stability of Lyophobic Colboids (Elsevier, Amsterdam, 1948). [21] H.C. Hamaker, Physica 4 (1937) 1058. [22] G. Kortism, Treatise on Electrochemistry (Elsevier, Amsterdam, 1985). [23] SW. Provencher, Computer Phys. Commun. 27 (1982) 213. [24] S.W. Provencher, Computer Phys. Commun. 27 (1982) 229. 125] P. Stépánek, in: Dynamic Light Scattering, The Method
[26] [27] [28] [29] [30] [31] [32] [33] [34]
[35] [36] [37] [38] [39] [40]
[41]
and some Applications, Ed. W. Brown (Oxford University Press, 1993) p. 175. J.A. Marqusee and J.M. Deutch, J. Chem. Phys. 73 (1981) 5396. D. Petsev and D.N. Denkov, J. Colloid Interf. Sci. 149(2) (1992) 329. D. Petsev, D.N. Denkov and K. Nagayama, Chem. Phys. 175 (1993) 265. D.N. Denkov and D. Petsev, Physica A 183 (1992) 462. B. Beresford-Smith, D.Y. Chan and D.J. Mitcell, J. Colbid Interf. Sci. 105 (1984) 216. S.B. Dubin, N.A. Clark and G.B. Benedek, J. Chem. Phys. 54 (1971) 158. V. Mikol, E. Hirsch and R. Giegé, J. Mol. Biol. 213 (1990) 187. A.J. Sophianopoubos, C.K. Rhodes, D.N. Holcomb and K.E. van Holde J. Biol. Chem. 237 (1962) 1107. W.H. Press, B.P. Flannery, S. Teukolsky and W.T. Wetterling, Numerical Recipes (Cambridge University Press, 1990). D.C. Phillips, Sci. Am. 5(215) (1966) 333. A. Watillon and AM. Joseph-Petit, Disc. Farraday Soc. 42 (1966) 143. R.H. Ottewill and J.N. Shaw, Disc. Faraday Soc. 42 (1966) 154. J. Visser, Adv. Colloid Interf. Sci. 3 (1972) 331. B.V. Derjaguin, Y.I. Rabinovich and NV. Crusaev, Nature (London) 272 (1978) 313. M.D. Soumpasis, in: Computations of Biomolecular Structures, Eds. M.D. Soumpasis and T. Jovin (Springer, Berlin, 1993) p. 223. M.D. Soumpasis and Y. Georgalis, in preparation.