1934
Langmuir 1994,10, 1934-1947
Calculation of Surface Polarity Profiles from the Experimental Hydration Force Data Stephan Kirchner and Gregor Cevc* Medizinische Biophysik, Technische Universitat Miinchen, Klinikum r.d.I., Ismaningerstrasse 22,D-8000Miinchen 80, E. U. Received July 30, 1993. In Final Form: March 11, 1994”
A method for the evaluation of the interfacial polarity profiles is introduced based on a previously proposed model (ref 10). This numerical procedure depends on the knowledge of a few, reasonably wellknown bulk solvent properties-such as the static and high-frequency dielectric constant. Water correlation decay length is argued to affect the results of the data analysis only little as long as its value is sufficiently smaller than the effective width of the interfacial region. The corresponding analysis of the hydration force data for the phospholipid bilayers shows that the interfacial polarity profile for such systems can be approximated well by a simple exponential function. The latter is uniquely defined in the region covered by the experimental data. However, without additional assumptions no extrapolations in the experimentally unexplored region are possible. In order to learn more about the molecular mechanism and the range of the surface hydration, as well as about the hydration-induced changes in the interfacial structure, more and more detailed experimental hydration force data should be collected.
Introduction Hydration plays an important role in many systems, ranging from the clay suspensionsto the biologicaltissues. In particular, the so-calledhydration force has been shown to govern the functional and packing properties of lipid multibilayers,’ polynucleic2or polycarbohydratestrands? and even proteins? In spite of this, the detailed mechanism of molecular hydration is only partly The same is true for the range of interfacial hydration which has been reported to vary between 0.16 and 0.45 nm,7 or even more.8 The first reason for this is the lack of detailed and practically useful models for the description of interfacig hydration. The other is probably the intimate coupling between the interfacial and the solvent structure; this is alsothe reason why the surface packing constraints (“steric phenomena”) may affect the range and the strength of so-called “hydration force”. We have recently shown how the effects of finite interfacial thickness can be incorporated into a meanfield theory of hydrati~n.~JOThe necessity for such a correction was clearly demonstrated in recent publication
by Marrink and colleagues,11who have also confirmed the basic correctness of our theoretical approach. All these studies reveal very strong effects of the interfacial structure on the structure and properties of surface-associatedwater. The distribution of surface polar residues, for example, seems to dramatically effect and sometimes dominate the spatial profile of hydration force near a hydrophilic surface. This effect has been suggested by uslo and others12to be particularly strong when the interfacial thickness exceeds appreciably the bulk solvent correlations length. In ref 10 a rationale has also been introduced for the evaluation of the interfacial hydration potential $&)-and thus of the hydration pressure-by relating these functions to the interfacial polarity profile, pP(x),l3 through
Jm
p,(x’)
exp[(x - x’)/Al dx’) (1)
Abstract published in Advance ACS Abstracts, May 15, 1994. (1) Rand, R. P.; Parsegian,V. A. Hydration forces between lipid bilayers. Biochim. Biophys. Acta 1990,988, 351-377. (2) Rau, D. C.; Parsegian, V. A. Direct measurement of temperaturedependent solvation forces between DNA double helices. Biophys. J. 1992,61, 260-271. (3) Rau, D. C.; Parsegian, V. A. Direct measurement of forces between linear polysaccharides xanthan and schizophyllan. Science 1990, 249,
A is the decay length of the (assumedly exponential) intrasolvent (this is, solvent structure) correlations. eo is the permittivity of free space. em (1.8I tmI 4.9)14is the short-wavelength dielectric constant. If the magnitude and the distribution of all local excess charges on the individual polar segments of the molecules are known, the corresponding polarity profile can be calculated directly. For this purpose, a local excess charge density for each polar group or atom, r, is first evaluated from the corresponding entity charge, Q,,and volume, V,.
(4) Kornblatt, J.,Hoa, G. H. B.; Heremans, K. Pressure-induced effects on cytochrome oxidase: the aerobic steady state. Biochemistry 1988,27, 5122-5 128. (5) Westhof, E., Ed. Hydration of Biological Macromolecules; Macmillan Press: New York, 1993. (6) Cevc, G.; Hirth, M. Spatial Hydration Profile of the Organic Materials at Submolecular Resolution. Submitted for publication. (7) Lyle, I. G.; Tiddy, G. J. T. Hydration forces between surfactant bilayers: An equilibrium binding description. Chem. Phys. Lett. 1986, 124,432-436. (8) Israelachvili, J. N. Intermolecular and Surface Forces, 2nd ed.; Academic Press: London, 1992. (9) Cevc, G. The molecular mechanism of interaction between monovalent ions and polar surfaces, such as lipid bilayer membranes. Chem. Phys. Lett. 1990,170,283-288. (IO) Cevc, G. Hydration Force Depends on Interfacial Structure of the Polar Surface. J. Chem. SOC.,Faraday Tram. 2 1991, 87, 2733-2739.
(11) Marrink, S.-J.; Berkowitz, M.; Berendsen, H. J. C. Molecular dynamics simulation of a membrane-membrane interface: The ordering of water and ita relation to the hydration force. Langmuir 1993,9,31223131. (12) Venable, R. M.; Zhang, Y.; Hardy, B. J.; Pastor, R. W. Molecular dynamics simulations of a lipid bilayer and of hexadecane: An investigation of membrane fluidity. Science 1993,262, 223-226. (13) Polarity profile, in the most general sense, corresponds to the spatial distribution of water binding sites (weighed so as to account for their different affinities and accessibilitiea)in the direction perpendicular to the surface. In the more specialized model used in this work, this function corresponds to the distribution along the surface normal of the local excess charges. (14) Komyshev, A. A.; Leikin, S. Fluctuation theory of hydration forces: The dramatic effectsof inhomogeneousboundary conditions.Phys. Reu. 1989, A40,6431-6438.
1278-1281.
0743-7463/94/2410-1934$04.50/00 1994 American Chemical Society
Calculation of Interfacial Polarity Profiles The standard pp,r = Qr/ V,relation is used for this purpose. Individuallocal excesscharge values can stem from suitable computer simulations or, even better, from direct measurements.16 (If required and known, the differential physical accessibilityto water, 0 Ictr I 1,of the individual polar groups can also be accounted for: ppJ = ctrQr/Vr.) The required volume values are deduced from Pauling radii after correction for the effects of thermal smearing. Finally, the integral surface polarity profile is determined by summing up all contributions from the polar residues and by the subsequent surface averaging
The dash signifiesthat only such polar groups are included in the above sum which are not compensated directly by mutual intrasurface charge-transfer processes. The latter may arise from intrasurface hydrogen bonds or from other direct intrasurface interactions. (This means that only polar residues free to bind water are important for the surface hydration.) Equation 1is based on the assumption that hydration potential of any polar surface is the primary source of interfacial hydration. It presumes, moreover, that this potential can be described within the framework of a onedimensional mean-field, nonlocal electrostatic model in which the value of A is a separation and position independent quantity. Equation 5 from ref 10provides a basis for the calculation of hydration force profile from the known interfacial polarity profile pP(x).l6 To date, however, direct information about the latter is not available. In order to calculate hydration force as a function of the interfacial thickness, ad hoc assumptions about the distribution and the magnitude of surface polarity sources are thus as yet necessary. In ref 10,for example, such sources were taken to stem from the density of the surface local excess charges. The latter were taken to decay exponentially with separation from the polar-apolar boundary. Alternatively, the surface polarity function was assumed to be of Gaussian type.17 In this contribution we try to resolve the ambiguity of interfacial polarity functions. For this purpose, we first reverse the procedure described in previous paragraph. Next, we show how the form of interfacial polarity profile can be deduced from the spatial dependence of "hydration force" data.1° Finally, we reconstruct the interfacial polarity profiles of several lipid membranes in contact with water by adjusting the parameters of our model so as to ensure optimal agreement between the appropriate fitting functions and the experimental data. Last but not least, we discuss which functions should be used for reliable and fast analyses of the hydration force data. Our new method requires no a priori knowledge about the function pp(x) or of the intersolvent correlationslength, A. This method thus minimizes the danger of a biased or prejudiced data analysis. Function pp(x) can be written as a series expansion in terms of some full set of (15) Kirchner,S.; Cevc, G. Hydrationof Polar Interfaces, a Generalised Mean-Field Model. J. Chem. SOC.,Faraday Trans. 2, in press. (16) In reality,the interfacial polarity profile is a function of all three spatial variables,p p ( x , y p ) . This latterfunction,to a good approximation, may be replaced by ita one-dimensional area average, however, when the interfacial thickness is greater than water correlations decay 1ength.M For the phospholipid bilayers this requirement appears to be fulfilled. (17) Equation 6 in ref 10 is partly in error. The first and the third exponent in brackete (the "pure"terms) lack a factor of 2 in front of x , owingtothe factthat in theexpreasionsfor the hydrationpreasure2p(x=d/ 2) should be used, to account for the membrane interaction in the superposition approximation.
Langmuir, Vol. 10, No. 6, 1994 1935 orthonormal functions. Eigenfunctions of a self-adjoint operator, that emerges in the mean-field theory of hydration, provide a suitable expansion basis for this. They ensure a rapidly converging data analysisand also simplify numerical calculations. To carry out this conceptualprogram, it is first necessary to reformulate the theory of hydration as a boundaryvalue problem, however. Electrostatic Mean-Field Theory of Hydration In refs 10 and 15, the hydration potential of an interface with a finite thickness d , was calculated within the framework of a generalized Poisson-Boltzmann model of the interfacial hydration. The corresponding (Helmholtz) free energy was then determined by an imaginary charging process of the interface according to standard procedures.lam Finally, the interfacialhydration pressure was derived from the thermodynamic relation p h = -(l/ A)aFh/ad,, where A denotes the total surface area. The justification of such a somewhat cumbersome calculation is given elsewhere.I6 It is based on the formal generalization of the original MarEelja-RadiE-Gruen model of surface hydrati0n.'~*~~*~3 The latter is given a different meaning, however: hydration in our model is not a manifestation of simple dipolar polarization (macroscopic electrostatics) but rather of direct water binding or the formation of water chlatrates around the surface polar residues (quantum-mechanical electrostatics).~~~~ Surface-induced water perturbation, consequently, is treated as a scalar rather than as a vectorial quantity. A more complete theoretical discussion is given in ref 15. In the linear approximation this is done rather directly by using Landau or Poisson-Boltzmann model and by introducing an extra term in the free energy functional. This then allows for the surface hydration potential (or field) For the model studied in this work, this corresponds to introducing a term that is a function of the spatial distribution of structural, membrane-bound charges pp(x)9J0921 into the free energy functional. The free energy of asingle hydrophilic solvated interface within the framework of such model is then written as
where the role of order parameter is now played by the hydration potential & and parameter K is defined as
K=
A2
coem(l
-4 4
(3)
e is the static dielectric constant. pp(K) is the interfacial polarity profile. In our electrostatic model, the latter corresponds to the distribution of local excess charges on each hydrophilic surface.26 (18) Gruen, D. W. R.; Marblja, S. Spatially varying polarization in Faraday Trans. 2 1983, 79,225-242. water. J. Chem. SOC., (19) Cevc, G. Molecular force theory of solvation of the polar solutes. The mean-field solvationmodel,ita implications and exampleafrom lipid/ water mixtures. Chem. Scr. 1985,25, 97-107. (20) See ref 15. (21) Kirchner, S.; Cevc, G. Derivation of a Generalized Surface Hydration Model, a Generalised Mean-Field Model. J. Chem. SOC., Faraday Trans., in press. (22) MarEelja, S.; RadiE, N. Repulsion of interfaces due to boundary water. Chem. Phys. Lett. 1976,42, 129-132. (23) MarEelja, S. Structural contribution to solutesolute interaction. Croat Chem. Acta 1977,49,347-358. (24) Cevc, G.; Podgornik, R.; &ki, B. The free energy, enthalpy, and entropy of phospholipid bilayer hydration and their dependence on the interfacial separation. Chem. Phys. Lett. 1982, I , 193-197.
1936 Langmuir, Vol. 10, No. 6, 1994
Kirchner and Ceuc
In this approach, the free energy per an isolated, partly hydrated membrane is given by
d$2 is the thickness of the (absorbed) water layer. By simple differentiationof eq 2 the interfacialhydration pressure of an isolated membrane, hydrated only between zero and the distance d , is found to be
z$h(d)&p(d)l (5) In the thermal equilibrium, the spatial profile of the hydration potentialof a fully hydrated surfacemust satisfy the Euler-Lagrange equation correspondingto eq 2. One can thus write with boundary conditions
$h(x==) = 0 (7) Solutions to these equations thus define the actual form of the interfacial hydration potential of an isolated, noninteracting membrane. For an arbitrary polarity profile this solution can be deduced formally from the corresponding Green's functions 0
$h(x)
= -JG(x,y)&p(y)
dy
(8)
0
Spatial dependence of the hydration pressure near a hydrophilic surface in such an approximation is thus strongly influenced by the surface polarity profile. To estimate the repulsion between two interacting membranes, a linear-superposition-approximationcan be used.26 Disjoining pressure in such an approximation is thence identified with the pressure difference between infinity and dw/2of two opposing "single membranes" (9)
This gives values which are probably somewhat higher than the correct result, owing to the separation dependence of the surface hydration potential. Based on the corresponding results for two surfaces interacting via Coulombic force26we expect, however, that eq 9 is in error by a factor of 2 or less. (We refrain from presenting exact results for the interacting membranes. The reason for this is that the separation dependence of the surface hydration potential is probably much smaller than the corresponding variation of the interfacial polarity profile, which is as yet unknown.) The geometry of our present model is illustrated schematically in Figure 1. In panel A, the case of an isolated hydrated polar surfaceis shown; in the lower panel, B, the situation of two interacting surfaces in a box of finite width b is shown (superpositionapproximation, see further discussion). The surface polarity profile, pp(x), and the interfacial hydration potential, h ( x ) , in the linear approximation (25) It can be shown relatively easily21 that 4 really describes the distribution of the surface-associated local excess charges. Constant K then takes into accountthat hydration potential, h,in the approximation used in this work, is connected with this charge distribution through a generalized, non-local-electrostaticPoisson-Boltzmann equation. (26) Verwey, E. J. W.; Overbeek, J. T. G. Theory of the Stability of Lyophobic Colloids;Elsevier: New York, 1948.
x = o
x =d42
x = b
1:
Figure 1. (A) Schematic representation of the model used in this work. Box width, b, is always at least slightly (in our calculationsbyO.1 nm) greater than the "half-separation" between the surfaces a t equilibrium, dw/2. The shaded area shows the contribution of to the free energy of hydration used for the evaluation of the hydration pressure in the linear superposition approximation. (B) In the linear superposition approach hydration potentials from both surfaces are assumed to be independent of the interfacial separation d. Hydration pressure in this approximation is then proportional to the magnitude of the "excess potential", 'h shown by the shaded area.
are both defined in the interval x = 0 to x = b, the latter value correspondingto infinite separations. (The reasons for introducing parameter b are discussed in more detail in Appendix A.) For the calculation of the hydrationdependent repulsion the hydration potential difference between the midplane (dw/2)and zero is considered (Figure 1). Implicit in our approximation is the overlap of the polarity profiles of the interacting polar surfaces for dw/2 smaller than the spatial extension of pp. This overlap can be kept small, however, if b 0. Even with an overlap possible, however, the hydration potentials used for the evaluation of "hydration pressure" do not interfere directly.n This means that mechanical contributions are excluded from our calculations: surface polarity profile is considered to be independent of the membrane distance d,.= From the experimentalist's point of view, the most interesting applicationof the equations introduced in this work is for the determination of the interfacial polarity profiles from the experimentalhydration-forcedata. Such a determination is a difficult task, however. The "testfunctions" for the analysis of polarity profiles should namely cause no bias to a special class of functions. In our analytical approach, pp(x), consequently, is taken to be a
-
(27) We believe that the partial overlap between the 'interfacial regions",and thus some overlap of the corresponding polarity profiles, is possible. After all, the interfacial polarity profiles are averaged probability distributions. Their overlap therefore does not mean that two lipid heads will occupy the same volume at the same time. Rather than this, it implies that two lipid headgroupswill partly share the same accessible region. Similar is the situation, for example,with lipid chains from both bilayerhalvessharing part of the membrane interior into which they penetrate. Our picture differs, however, from the one advocated by Israelachviliand Wennerstr6m.B Whereas these authors consider direct forces of the steric origin, we only deal in this work with hydrationmediated forces between the interacting headgroups. (28) To include the mechanical (entropical)variabilityof Pp(x,d-), one should abolish the simple model proposed in this paper. Instead, one should calculate the partition sum of a system that includes both the mechanical as well as hydrationeffects. In sucha description,for example, the ansatz of Israelachvilliand Wennerstrem can be combined with the mean-field theory of hydration by replacing the hard-core interaction of the former model by a hydration-dependent Hamiltonian. The resulting partition sum could formally look like: Z = Sod., dZ1 J0d.r" d 2 2 exp{-(l/ kT)ffhydm(Zlp2,dw)lwith ffhydm = -rK J'b(rrls2,dw)~,(rp192,d2) dr, h,and pp being connected through eq 8.
Langmuir, Vol. 10, No.6, 1994 1937
Calculation of Interfacial Polarity Profiles sum of such basis-functions that can represent any arbitrary function. The easiest way to construct an initial guess for the function p&) is thus to use a box model with a constant box width and a variable box height. In this case, the interfacial hydration pressure can be calculated analytically. Such an "optimization" does not lead to reasonable fits, however. This is chiefly owing to the fact that the interfacial pressure decays exponentially with the solvent correlations length over the width of each box. The resulting model parameters in such an approximation, consequently, are highly variable and strongly dependent on the starting model parameters (not shown). To improve on this, one can work with a set of orthonormal, such as Laguerre or Hermite, polynomials. Unfortunately, this approach does not perform well either. This is chiefly due to the numerical reasons. During each optimization run, repeated numerical integrations and subsequent numerical differentiations are done. This introduces many errors. We did not succeed in obtaining reasonable polarity profiles by fitting the measured hydration force data in this way. Solution of the Fit-Problem by Diagonalization of a n Approximated Green's Function A good way to circumvent the mentioned difficulties is to diagonalizeGreen's operator by the procedure described in Appendix A. This having been done, any external hydration potential can then be expressed as a series expansion in terms of the eigenfunctions C, which fulfill the boundary conditions given by eq 21.% Starting with the polarity profile m
one gets for the interfacial hydration potential m
and finally for the hydration pressure
m
where*
Cn(d/2) = (2/bI1/' sin[nad/(2b)]
(14)
(29) Correspondingdiscusion and evidence are given in every textbook on ordinary differential equations that also includes a chapter on the boundary-value and eigenvalue problems; the work cited (ref 30) is j u t one example for this. From a physicistspoint of view, this resembles the formalism of quantum mechanica, transformed into another language. Actually, the idea for our present approach waa guided by this analogy. (30) Walter, W. Gew6hnliche Differentialgleichungen;Heidelberger Taachenbacher 110; Springer, Berlin, 1976. (31)As explainedin AppendixA, there is some freedom of choice when and fo. Those given here were used only for choosing functions convenience.
The calculation of derivativesCn',Io', and fd, that appear in eq 12, is trivial. In terms of these functions hydration pressure can be expressed as a function of hydration potential that fulfills the following boundary conditions: $(x=O) = '15'0 and $(x=b) = $(x=m) = 0. Quantities AI, A2, ...,QOand p then remain unknown. (The appearance of b in eqs 9 through to 13 is caused by the fact that an infinite system must be replaced by a finite system of nearly arbitrary, constant width b; without such a modification no sum expansion is possible. This is tantamount to replacing the exact Green's function with an approximate one. Such an artificial restriction does not influence the final model results, however, as is shown in the results section, as long as the box-size b is much greater than the characteristic hydration decay length A.) Obviously, b must be chosen so as to be greater than the maximum value of dw/2. The position of x = 0 plane can then be located anywhere before the first data point (see further discussion). Expression 12 only contains simple sums and products. Its value, consequently, can be calculated rapidly and easily. All integrations and differentiations, which are relatively difficultto handle numerically, in this expression are replaced by the far more convenient summation terms. Infinite series in 12, moreover, can be approximated by a finite sum without a significant loss of accuracy. This is owing to the fact that our expansion basis is optimally suited for the problem tackled. The unknown Coefficients AI, ...,AN, '15'0,p are determined easily by a nonlinear leastsquares fit to the experimental data. A how-to-dosummary of the entire procedure is given in Appendix B. A big advantage of the proposed method is that it does not depend on a priori information about the actual form of the interfacial polarity profile. The reason for this is that any continuous function in the interfacial region, this is, between zero and the outer boundary 6 , can be approximated by means of a Fourier expansion with nearly arbitrary accuracy. Every such function, consequently, can be parameterized in the form of Fourier-expansion coefficients. Mathematical relation between the interfacial polarity and hydration pressure profiles is defined uniquely for each value of the solvent correlations length A (cf. Appendix A). It is clear, however,that the proposed calculation scheme is only an approximation that contains an inherent inconsistency. The free energy of a partly hydrated system is namely calculated by means of functions pp and hthat are defined for a fully hydrated system. One consequence of this is that the electroneutrality condition is only satisfied for dw/2 = b, since only then the local excess charges on the water molecules, -JodJ2 $h(X) dx, precisely cancel the local excess charges located on the surface polar residues, -Jodw/2 p p ( x ) dx. We further comment on this in thediscussion but let it be said already that this difficulty only leads to small errors for the reasonably hydrated systems.
Kirchner and Cevc
1938 Langmuir, Vol. 10,No. 6,1994 1
-50
' 0
0.2
0.4
5
0.6
0.8
'
[nml
-40
'
0
0.1
0.2
0.3
0.4
0.5
0.5
0.1
0.8
0.9
x [nm]
Figure 3. Effect of the changing number of the adjustable model parameters on the quality of the data analysis. Surfacehydration potential was assumed to have a Gaussian form, shown in Figure 2C. Deviation between the original and fitted profiles is shown for the different number of model parameters: (dotted)9 parameters; (dashed) 10 parameters; (dashed-and-dotted)11 parameters; (full) 12 parameters.
x [nm]
Figure 2. Effective surface polarity profile (A and C)and the corresponding disjoining pressure profiles (B and D) calculated within the framework of a generalized electrostatic model of surfacehydration. The starting data (shown as +) were expressed in terms of 64 Fourier coefficients; for the final numerical data analysisthis number was reduced to 10 (8expansion coefficients, p, and PO). The starting surface polarity profile, represented by crosses,was reconstructed by the described analytical procedure nearly perfectly. The difference between the results of the fit ('Gaussian" curves) and the starting data (+) is always quite small (see the low-amplitude curves). (A, C) bo = 1000 e~p-0~6~(d.r0.2)/0.161~; (B,D) bo= 1()00exp-0.6[(d.r0.4)/0.1l1.
Results In order to demonstrate the power and usefulness of our analytical approach, we have first tried to determine the interfacial hydration potentials from a set of the simulated hydration force data. We have thus calculated, via eq 12, the hydration pressure of two surfaces with a Gaussian distribution profile. This pressure was then fitted by using similar equations with a much smaller number of characteristic parameters A,. The calculated surface polarity profile, determined by the optimal parameter values, was finally compared with the original interfacial polarity profile. The results of these operations are illustrated in Figure 2.32 They vindicate the conclusion that the described analytical procedure ensures a very good agreement between the starting and deduced surface polarity profiles. Small deviations in picture 2c get even smaller if only one additional adjustable parameter is allowed for. We do not show these more accurate results in order to facilitate a direct assessment of the reliability of our analytical procedure for the model and real surfaces. (32) To calculate the pressure from a given profile k = f&f), we have first calculated the Fourier coefficients of f ~ ( d -) A&l) with Aw(d) = fcC0). sinh(?;(b- d)!/sinh(2b) by numerical integration (i.e. p = 1). To obtain the interfacial hydration pressure, eq 12 was then used.
Deviations in picture 2c originate from the fact that the sine function of the highest order chosen (here nmax= 8) is still too coarse for a good reconstruction of the sharp tip of the rather narrow starting Gaussian function. For a good profile reconstruction a reasonably high number of the adjustable model parameters thus should be chosen. This is documented for the narrower of the two illustrated Gaussian curves in Figure 3. Another aspect of the finiteness of our system deserves mentioning as well. If only a fraction of the water-filled region is known experimentally, our analysis can only provide good results for such an experimentally studied region. For the same reason, the precise value of the boxwidth b does not affect the accuracy of the experimental data analysis. Figure 4 demonstrates this clearly. The reconstructed surface polarity profile corresponding to one of the Gaussian curves shown in Figure 2 is virtually the same for b = 1,1.2, and 1.5 nm. An extra peak is observed for b = 1.5 nm only. This is caused by the fact that no data are available beyond 1 nm. The reconstructed polarity profiles thus lack constraints in this region and, consequently, are unreliable. This immediately shows one limitation of the analytical method proposed in this work: our approach can only be used for analyzing the experimental data and not for extrapolations. When data are lacking in some region, the reconstructed surface polarity profile in such a region must be fixed through ad hoc assumptions or else remain obscure. Most figures presented in this work corroborate this fact. After having tested the potential limitations of our new procedure for the analysis of hydration force data, we have also calculated the polarity profiles of several real surfaces by the same method. The published hydration force data for the dipalmitoylphosphatidylcholine (DPPC) multilamellae at 25 and 55 "C were selected for this purpose, being representative of the lipid bilayers in the gel- (L'o) or fluid-lamellar (La)phase, respectively. Near equilibrium, the net force between lipid bilayers falls steeply to zero. The reason for this is van der Waals
Langmuir, Vol. 10, No. 6, 1994 1939
Calculation of Interfacial Polarity Profiles 20
A
0
8.-
0
I 0
0.2
0.4
0.8
0.8
1
1.2
1.4
= [nml Figure 4. Effect of changing box width, b, on the results of the theoretical analysis of the hydration force data. Gaussian curve shown represented in Figure 2A was used as the starting, “experimental”data set. Box width was then set to b = 1,1.2, and 1.5 nm, and the fitting procedure was done with very similar results by using 10 adjustable parameters. attractions which becomes progressively more important at large separations and finally fully compensates the interfacial repulsion. For a complete description of the interfacial interactions one should, therefore, consider hydration as well as van der Waals force, at least. Fluctuation forces for the fluid-phase bilayers could also play some role. In the d, 0 limit the contribution from steric forces, in principle, could play some Such forces can only manifest themselves in “pure” form in the absence of headgroup solvation. As soon as they start to depend on the water-headgroup interactions, they become a sort of interface-dependent “hydration force” and are, as such implicitly included in our model. In this work we are primarily interested in the surface polarity profiles. These are affected less by the van der Waals attraction than by the solvent-mediated force, especially at distances smaller than the equilibrium interfacial separation. In order to avoid dealing with dispersion forces, we have therefore excluded all data pertaining to the largest observed separations from the data analysis. Only data points falling on the curve with a characteristic quasi-exponential decay were therefore analyzed.34 Equation 12 contains numerous unknown parameters. In order to be applicable to the rather sparse real data, we have thus first constructed sets of the corresponding but more numerous quasi-data. This was achieved by splinefitting the original hydration-force d a h 1 Equation 12 was then used to analyze the resulting spline curves. Spline-fitted, quasi-experimental but “higher quality”, data are shown as dots in Figure 5. The results of the corresponding data analyses are represented as solid lines. It is clear that the two sets of data are nearly identical.
-
(33) McIntosh, T. J.; Magid, A. D.; Simon, S. A. Cholesterol modifies the short range repulsive interactions between phosphatidylcholine membranes. Biochemistry 1989,28,17-25. (34) If van der Waals force is allowed for, by using a Hamaker constant value as measured in direct force experiments (H= (6-12) X lo-*’J) the analysis, based on the assumedly exponential polarity profile, suggests that the repulsion at large separations is qwi-exponential with a very short decay length of approximately0.1nm. Model results for the shorter separations are only modified in a trivial manner, however.
h
9
.2ti 4
04
1.8
2.2
2.8
3
4
d , [nml
Figure 5. Measured (symbols)and spline-fitted (curve) hydration pressure profiles for the gel-phase dipalmitoylphosphatidylcholine multi-bilayers at 25 O C (A) and for the similar multilamellar stacks of DPPC in the fluid phase at 55 O C (B). (Data stem from ref 1.) In all simulations, eight adjustable parameters were used (6 expansion coefficients, p, and 90). Solvent correlations length value was taken to be A = 0.075 nm.
(Allanalyses were based on functions A d and l o defined in previous section. This is only a matter of convenience, however. Other forms of these functions give similar results.) To determine the optimal model parameter values by the nonlinear data regression we have used the Marquard algorithm.% The bulk intersolvent correlationslength was taken to be A = 0.075 nm? unless stated otherwise. This value gives the apparent decay length of the hydration phenomena in the solutions of small polar organic molecules for which no interfacial or steric effects are possible. We thus prefer to use this estimate rather than previous less directly measured values (see also further discussion). It is noteworthy, however, that even using more conservative estimates such as A = 0.15 nm or even 0.18 nml leads to qualitatively similar results. Box width was normally taken to be marginally greater than half of the equilibrium distance between the interacting polar surfaces ( b = d d 2 + 0.1 nm), but smallvariations in this value did not affect the results of our model appreciably.% Parameters and p were permitted to vary during each fitting procedure. The latter parameter, which has no real physical meaning, could also be held constant. Fitting procedure would then rely just on the expansion coefficients A,,. However, the use of a variable parameter p offers the opportunity of fitting an approximately
*,-,
(35) Bevington, P. R. Data Reduction and Error Analysis for the Physical Sciences; McGraw-Hik New York, 1969. (36) In light of the boundary conditions 8, it is better to w e a box of width b i:dw+,/2. In reality, of come, b = dw/2,owing to the van der Waals attraction between the lipid membranes.
Kirchner and Cevc
1940 Langmuir, Vol. 10, No. 6,1994 14000
T
,
I\ Iy
i
A
I /
I
I 0
0.2
0.6
0.4
0.E
8
1
x [nin]
Figure 6. Effective surface polarity profile of dipalmitoylphosphatidylcholine bilayers at 25 O C as deduced from the experimental data shown in Figure 5A by means of the theoretical
analysis described in the main text. Various surface polarity profiles correspond to the same experimentaldata set. Different curves were obtained by choosing various startingvalues for the adjustable model parameters. Quasi-exponentialfit (which is most linear in the logarithmic representation) has the lowest +value and is thus closest to the global optimum. Region for which the experimental data are available is enclosed by two vertical lines and the corresponding polarity profiles are given as full curves. Dashed lines correspond to the ill-definedregions.
exponentially decaying experimental data with a quasiexponential function. This leaves the task of reproducing the differences between the quasi-exponentialand the real surface polarity profiles to parameters An solely. If this difference is small, An values are also small. By keeping p # 0 the effective number of the adjustable model parameters is thus minimized, provided that the experimental curve contains no fine, “high-frequency” features. Typical results of the corresponding fitting procedure are shown in Figures 6 and 7. They illustrate the consequences of starting-parameter changes. They also suggest that various calculated polarity profiles may differ dramatically in the experimentally unexplored vicinity of the hydrophilic surface if unreasonable surface values are chosen. This notwithstanding, all fits are essentially identical in the experimentally secure region covered by the measurements. The agreement between measured and calculated “hydration force”data in this range is essentially perfect (see the minimal deviation between the curves and symbols in Figure 5). Final results in Figures 6 and 7, which pertain to the dipalmitoylphosphatidylcholine multi-bilayers, thus do not represent an incidentally reached local minimum but a true global minimum (cf. Appendix A). From the logarithmic plot of the calculated surface polarity profiles, one can see that in the region that is reliably backed-up by the experimental data the function
Figure 7. Surface polarity profile of the dipalmitoylphosphatidylcholine multi-bilayers at 55 O C reconstructed from the corresponding “hydration force” data as described in the text and in captions to 4 and 6. The specificvalues of the (arbitrary) model parameters were (from top to bottom) (I.C and \ko) as follows: (A) 0.211 and 13095,8.6 X lo7and-373,0.285 and 4119, and 0.034 and 32487; (B) 0.205 and 61010,0.214 and 52545,and 0.285 and 23015. pp is
nearly exponential, for the dipalmitoylphosphatidylcholine multi-bilayers at least. Exponential approximation introduced in our previous publication for the description of the interfacial polarity profileslo thus has been a good choice but a suitably chosen Gaussian value would yield similar results. This is certainly true not too close to the membrane surface. The calculated interfacial polarity profile depends on the precise value of interwater correlations length used. In order to estimate the effect of this model variable on the final model results, we have, therefore, reanalyzed the same experimentaldata by using severalpossible,tentative A values (Figure 8). The results of these calculations suggest that the water correlations length affects the final model results only little as long as its value is sufficiently smaller than the apparent decay length of hydration force. When this condition is fulfilled, various A values merely decrease the absolute height of the surface polarity profile slightly and, as is shown in the logarithmic plot, increase the slope of the p, vs d data a little. These results are not unexpected. In the exponential approximation’o the interfacial hydration pressure is found to consist of three different terms. Two of these decay exponentially on the length scale of A and d,, respectively. The third term is mixed and depends on both these characteristic lengths. The overall influence of the solvent correlations decay length on the surface polarity profiles is normally quite small, however. This is owing to the fact that interfacial thickness is typically much larger than this decay length, d, >> A. The range of polarity profiles, in particular, is nearly independent of the assumed values
Calculation of Interfacial Polarity Profiles
Langmuir, Vol. 10, No. 6,1994 1941 1.4
T s 8
1
0.6 0.21 0
0.4
0.8
1.2
x
0
0’2
0.4
0’6
x
2
1.6
[nml
i
0’8
[nml
12
B
II
I 0
02
0 4
06 2
O B
0
08
[nml
Figure 8. Influence of the intrinsic decay length (i.e. ‘solvent correlationslength”),A, on the model-derivedsurface hydration potential and on the calculated surface polarity profiles of dipalmitoylphosphatidylcholine multi-bilayers. The original hydration force data were the same as in Figure 3. They stem from the measurements with DPPC at 25 O C . Eight adjustable parameters were used to get the shown fit (seealso texts to other figures). Box width was b = 1 nm throughout. Decay length values were A = 0.005,0.1,0.15,0.18,and 0.2 nm, from top to bottom, respectively. (Themeasured apparent ’hydration-force” decay-lengthvalue for this system is b d m a= 0.19 nm.9 The results of data analysis are very similar except for A = 0.2 nm. for the interwater correlations. In light of this interpretation we consider the solvent dependence of the experimental A values, reported by McIntosh and c0lleagues,3~ to be the result of the different swelling and/or mobility characteristics of the interfacial region in different solvents. (This is tantamount of assuming that d, = d,(solvent) in our present model.) In the hypothetical case in which the relevant intersolvent correlations length is much greater than the effective decay length of experimental “hydration force”, the polarity profile would have the form of a &function. Hydration pressure would then be determined fully by the former decay length value. For the dipalmitoylphosphatidylcholine bilayers in the gel phase this would require the range of correlations between the interfacially bound water to be greater than the experimentally determined value of 0.19 nm, however. In our opinion, this is rather improbable6 (see also further discussion). Hitherto, we have only considered different possibilities for the reconstruction of the surface polarity profiles based on the measured “hydration-pressure” data. It is equally possible, however, to calculate the interfacial hydration potential and hydration pressure from the given interfacial charge-distribution functions. (37) McIntosh, T. J.; Magid, A. D.; Simon, S.A. Range of the solvation pressure between lipid membranes: dependence on the packing density of solvent molecules. Biochemistry 1989,28, 7904-7912.
\
\ !
1
16
2.4
32
4
d, [nm]
Figure 9. Molecular structure of a typical phospholipid, phosphatidylethanolamine(top,atomic coordinatesfrom ref 38) and the correspondingexperiment-derivedspatial profiles of the interfacialpolarity (pp(x), second from the top),surfacehydration second from the bottom), and interfacial potential (ha(%), bottom). Atomic distribution in the hydration potential @-, top panel (open ellipses, H; filled ellipses, C; dashed ellipses, 0, circle, N; large circle, P) corresponds to the crystal structure data for dilauroylphosphatidylethanolamine.98 The corresponding experiment-derived intrfacial polarity profile (second from the top)was calculated as described in the text by combining the structural information with the local excess charge densities on the phosphorylamine crystals measured at 123 K m in order to estimate the excess charge densities on the glycerol moiety we have used the method described in ref 19. From such nominal charge density values the appropriate density profile was calculated by means of eq 18. Surface potential and surface hydration pressure profiles were obtained from pp(d)as described in Appendix B. One example for such a calculation is given in Figure 9. It shows the interfacial polarity profile obtained by the integration of the charge-distribution profile for the dilauroylphosphatidylethanolamine multibilayers (DLPE). This calculation relies on the crystallographic data input from ref 38 and on the use of the local excess charge values measured directly for the polar residues of the phosphoryl headgroup crystals by Craven and his colleague^.^^^^ Based on the results of our quantum-mechanical calculations made for the different temperatures with the commercially available computer simulations, we claim ~
~~~~
~
(38) Elder, M.; Hitchcock, P.; Mason, R.; Shipley,G. G. A refinement analysis of the crystallography of the phospholipid, 1,2-dilauroyl-DLphosphatidylethanolamine,and some remarks on lipid-lipid and lipidprotein interactions. R o c . R. SOC.London, A 1977,354, 157-170. (39) Swaminathan,S.;Craven,B. M. 1984. Electrostatic properties of phosphorylethanolamineat 123 K from crystal diffraction. Acta Crystallogr. 1984, B40, 511-518. (40) Weber, H.-P.;McMullan, R. K.; Swaminathan,5.; Craven,B. M. The structure and thermal motion of phosphorylethanolamine at 122 K from neutron diffraction. Acta Crystallogr. 1984, 840, 506-511.
Kirchner and Ceuc
1942 Langmuir, Vol. 10,No. 6,1994
that these charge values, which are the only directly measured ones reported to date, can also be used to illustrate the situation encountered at ambient temperatures. Local excess charges that we have used in this work to characterize the glycerol-backbone region are estimates deduced on the basis of structural ~imi1iarities.l~ The one-dimensional,low temperature surface polarity profile shown in Figure 9 was calculated by first projecting a single phosphatidylethanolamine molecule onto the membrane surface director. From the absolute values of local excess charges on the polar residues of such a molecule and from the spatial distribution of the corresponding segments, the net charge density profile along this director was then evaluated. This was achieved by summing-up all individual charge densities in two steps. First, each atom i in position ri was considered to be a homogeneously charged sphere. The radius of the latter was identified with the corresponding van der Waals radius Ri. The effective charge density of each group or an atom i was thus taken to be given by: p i = Qi/(47rR$/3). The calculated charge-density value was then assigned to the appropriate position on the membrane normal and divided by the lipid area to account for the lateral packing density effect. This gave a series of discrete charge values a t different positions along the bilayer normal. In order to allow for the thermal vibrations of individual groups, we have next smeared these individual values, in a Debye-Waller modus by using the measured width of the thermally-induced Gaussian distribution of the polar groups or atom positions (at 123 K), ai. Finally, all local excess charge value distributions were summed up to yield41
This procedure is very similar, for example, to the calculationof the net charge distribution in adiffuse double layer than stems from the presence of the individual ions next to a charged surface.
Discussion In the classical mean-field theory of charged surfaces in contact with an aqueous electrolyte, so-called diffuse double layer theory, the constraint of an infinitely narrow interface was first relaxed by Stern. This authofi2 has introduced the concept of an adsorbed ion layer that is also inherent in the Donnan-layer approach.43 Later on, a corresponding generalization of the Gouy-Chapman approach was proposed for the integral description of the charged surfaces with a thick interface.44 This work demonstrates how a similar approximation can be used to calculate the polarity profile of a hydrophilic surfaces in contact with an aqueous solution from the (41) Surface hydration potential shown in Figure 9C does not contain contributions from the surface bound water molecules. Such water molecules contribute to the s u m of the intrinsic, h,and polarization,$p, potentials,mwhich are determined with the surface ('dipole") potential measurements.3' It would be misleading, however, to include such contributions into h since the water-dependent potential part is a consequence rather than the cause of the surface hydration phenomena. The sum of such potentials is, in fact, essentiallynil, as hae been shown in the recent molecular-dynamics simulation of hydrated DPPC by Marrink and colleagues.11 (42) Stern, 0. Zur Theorie der elektrolytischen Doppelschicht. Z. Elektrochem. 1924,30,50&516. (43) Ohshima, H.; Ohki, S.Donnan potential and surface potential of a charged membrane. Biophys. J. 1986,47,673-678. (44) Cevc, G.; h k i , B.;Svetina,S.Theelectrostaticpotential of bilayer lipid membrane with a structural surface charge smeared perpendicular to the membrane-solutioninterface. An extension of the Gouy-Chapman double layer theory. J. Phys. Chem. 1981,85, 1762-1767.
experimental hydration force data. Simpler hydration models which only deal with the infinitely narrow interfaces are described in refs 14, 18-20, 22, 45, and 46. The results of our analysis are, of course, influenced by the reality of the following assumptions: First, we have postulated that solvent correlations length is not a function of the interfacial separation. This is probably true for the very short correlations lengths and intramolecular interaction modes considered in this work. Second, the water correlations decay length A was taken to be position independent. This postulate may, but need not, be true. Ita effect on the outcome of our analysis is rather small, however, as long as A
