J. Phys. Chem. 1991, 95,6351-6360
C. Two-Pulseand Tbree-Pulse X-Band ESEM Spectra of G1. The two-pulse ESEM spectrum of the glass G1 at 4.2 K is presented in Figure 9, while the three-pulsed ESEM of G1 at 4.2 K appears in Figure IO. Analyses of these modulated echo patterns are based on the appropriate theoretical expressionsIO for the two-pulse ( V ( 7 ) )and three-pulse (V( r ) ) modulated echo intensities for an electron spin (S = with axial g interacting For the local geometry of weakly with a nuclear spin ( I = close-packed ionic spheres envisaged in our model, ESEM spectral simulations based on spherical averaging" and used in the recent literaturell*Mwould be ideally suited. Specifically, our simulation procedure compared the experimental ESEM spectra with the calculated ones to give "best-fit" values for the Cu2+ electron~ O - ~ distance, O ~ P ~ r, the number of 207Pbnuclei, N , involved in this superhyperfine interaction, and the estimated value of the isotropic part of this electron-nucleus interaction, ab To account for the overall magnetization decay by magnetic spin-spin interaction, the calculated ESEM spectra were multiplied by a decay function of the form exp(C, + CIT C2P C37'), which was fitted by a nonlinear least-squares technique to the experimental decay curve." The modulation observed for the two-phase echo decay curve (Figure 9) is faint, indicating that the superhyperfine couplings between a Cu2+electron and close-lying Pb2+ nuclei in the first coordination shell ( 1 3 A) are very weak, if they are present. More interesting results are observed in the three-pulse study of the glass system (Figure 10). Here the decay curve is seen to be modulated by the presence of 207Pbnuclei in the vicinity of Cu2+. The presence of modulation only in the three-pulse echo suggests that Pb2+ is at slightly farther (>3 A) positions, conceivably in the second coordination shell, as would be expected from our X-ray diffraction results. The two dashed lines shown in Figure 10 are the trial-and-error computer simulations, which closely match the experimental echo decay, both computed assuming aim= 0.3 MHz. Figure 10b is the echo decay curve simulated on the basis of three Pb nuclei ( N = 3) surrounding Cuz+at a distance of 3.6 A in our spherical
+
+
(33)Kevan, L.;Bowman, M. K.; Naryana, P. A.; Boeckman, R. K.; Yudanov, V. F.; Tsvetskov, Yu. D. J . Chem. Phys. 1975, 63, 409. (34)Kevan, L. Acc. Chem. Res. 1987, 20, 1.
6351
averaging model. One observes, however, that the simulated echo modulations are too strong, particularly in the early-time ( T = 1 ps) region of the ESEM envelope. Further, trial-and-error procedures were tested to improve the overall fitting, such as the inclusion of 35Cland 37Clnuclear quadrupole and hyperfine interaction terms from distant chlorines, but these led to a suppression of the overall 3-pulse decay. On the other hand, the fitting in the early regions of the decay improved distinctly if we assumed superhyperfine interactions from four Pb nuclei ( N = 4) at a distance of 3.8 A. This is an interesting result, considering that we have assigned an X-ray PDF value of 3.8 to Pb-Pb or P W u scattering distances. The ESEM results therefore further confirm our earlier picture of the second coordination shell in G1. Such types of distance assessment where a direct chemical bonding is not involved are rather difficult to make by other studies. (For example, EXAFS studies2 and neutron diffracti01-1'~primarily yield information regarding the first coordination shell.)
Conclusion By combining the results of our X-ray diffraction, EPR and ESEM studies we have derived useful structural information from a novel 43Pb0-56PbCl2-CuCIZ ternary glass (labeled GI). The X-ray PDF data enable us to propose a general structural model for the glass in which octahedral building units of Cu04C12are predominant. EPR spectra of G1 at X- and Q-band microwave frequencies have been fitted, by line shape analysis procedures, to appropriately distributed spin-Hamiltonian parameters, suggesting that the Cu2+bonding characteristics in the glass structure are distributed. Among the distributed geometries, Cu2+ in the elongated octahedral unit, Cu04C12,appears to make the major contribution, even though other Cu2+ geometries such as pseudotetrahedral may be present to a minor degree. Our ESEM results suggest that Cu2+is surrounded by four Pb2+in the second shell a t a distance of 3.8 A.
Acknowledgment. P.R. is thankful to Professor Larry Kevan for making available pulsed-EPR facilities at the University of Houston. (35)Wright, A. C.;Grimley, D. I.; Sinclair, R. N.; Rao, K. J. J . Phys. (Paris) 1985, CB, 305.
A Simulation Study of Flexible Zwitterionic Monolayers. Interlayer Interaction and Headgroup Conformatlon M. K. Granfeldt and S. J. Miklavic* Physical Chemistry 2, Chemical Center, Box 124, 221 00 Lund, Sweden (Received: January 3, 1991) A simple model for two opposing monolayers of flexible zwitterionic amphiphiles is studied by means of Monte Carlo simulation. The model allows for molecule fluctuations with a component normal to the average hydrocarbon/water interface arising from the partial diffusion of the hydrocarbon chains out of the hydrophobic core. For a range of model parameters we investigate two properties of the system: the headgroup conformation (orientation) and the separation dependence of the interlayer force. From a comparison with spectroscopic data we find that the average "near-parallel" orientation of the zwitterion dipole is as much a consequence of the perpendicular fluctuations as of inter- and intralayer electrostatic interactions. The out-of-plane degree of freedom extends the range of a repulsive entropic contribution to the total interlayer force. The total force, which can be repulsive up to separations of 30 A, is found to have an exponential dependence on surface separation. At larger separations attractive correlation contributions dominate, providing a limit to the 'swelling" of the layers in excess water. The model produces trends qualitatively consistent with force measurements in systems of various phospholipids.
Introduction For some years now very large but short-ranged repulsive forces have been found to exist between charged and electroneutral and their origin has been the subject of intense specuTo whom correspondence should be addressed.
lation and debate ever since. These repulsive forces, over and above ordinary DLVO tYPe forces,3 have an exponential dependence on ( I ) Rand, R. P.; Parsegian, V. A. Eiochim. Eiophys. Acta 1989,988,351. ( 2 ) Israelachvili, J. N. Intermolecular and Surface Forces; Academic
Press: New York, 1985.
0022-3654/91/2095-6351$02.50/0 0 1991 American Chemical Society
Granfeldt and Miklavic
6352 The Journal of Physical Chemistry, Vol. 95, No. 16, 1991 the separation, 6, P(6)= Po exp(-b/A), with A of the order of a few angstroms, and are generally only present at small surface separations (b I30 A). These two often-quoted features, found to be similar for all the measured forces, have formed the basis of speculation of a unique mechanism. Recently, Rand and Parsegian' have reviewed the data on these so-called hydration forces, found in zwitterionic and charged phospholipid systems. This review should be considered side by side with measurements of forces between charged mica surfaces43 and with the forces measured between nonionic surfactant layers6 since all these systems find repulsive forces of some description, above that of ordinary double-layer interactions. If a unique mechanism is at work, then it would seem reasonable to expect that, to within some margin, these forces would respond similarly to changes in external conditions. In fact, this is not the case. For example, an increase in temperature in the phospholipid system results in an increased repulsion, in both range (A) and magnitude (Po),while it reduces the repulsion in the case of pentaoxyethylene dodecyl ether coated mica surfaces6 and even introduces a repulsive minimum gradually turning into an attractive potential well with increasing temperature. To the authors' knowledge the temperature dependence in the case of bare mica has not been studied. One other aspect, the addition of excess monovalent salt, is found to have no significant effect on the repulsive forces between phospholipid bilayers while the hydration force measured between bare mica is extremely sensitive to both cation type and concentration. Specifically, it is absent although at sufficiently low concentrations of most salt types and absent for all measured concentrations of an acid.4 Several explanations have been suggested to account for the measured repulsive forces. Two of these are an image charge repulsion' and surface-induced water polari~ation.~-~ Unfortunately, in seeking an analytical exponential repulsion, the original analyses have fallen short of achieving sufficient conviction mostly as a result of uncertain approximations and assumptions.loJ' Furthermore, the ideally suited work of Pashley4 did not at all support the idea of a water polarization mechanism. After initial onset of the short-ran e forces, Pashley's estimated decay length varied from 3 to 12 depending on salt type, suggesting more of a (hydrated) solute mediated repulsion rather than solvent mediated. Ii' this is a true understanding of the events of this system, wha: can then be said in the phospholipid case? In the forerunner of this paper" we sought to test the suggestion that interactions with image charges may be responsible for the measured repulsions. In a Monte Carlo simulation the phospholipid headgroup was modeled as a zwitterion located close to an interface between a high dielectric medium (water) and a low dielectric medium (the hydrocarbon core). The zwitterion was allowed two-dimensional translational freedom on a plane a certain distance away from the discontinuity and also allowed to rotate about the negative ion position. The force between two opposing layers was calculated for different locations of the plane of negative ions. We found in most of these calculations that although the image charge contribution is repulsive, the strong electrostatic interlayer correlations lead to an overall attraction. An asymptotic analysis of this or a closely related system has also shown that even if image charge contributions were included, the total pressure at large separations remains attractive.12 While this analysis and
A
(3) Verwey, E. J. W.; Overbeek, J. Th. G. Theory of rhe Stability of Lyophobic Colloids; Elsevier: New York. 1948. (4) Pashley, R. M. J . Colloid Interface Sci. 1981, 83, 531. Miklavic, S. J.; Ninham, 9. W. J . Colloid Inrerface Sei. 1990, 134, 305. ( 5 ) Israelachvili, J. N.; Pashley, R. M. In Biophysics of Water; Franks, F., Ed.; Wiley: New York, 1982. (6) Claesson, P. M.; Kjellander, R.; Stenius, P.; Christensson, H. K. J . Chem. Soc.. Faraday Trans. 2 1986, 82, 2735. (7) J a m n , Be.; Wennentrom, H. J. Chem. Soc., Faraday Trans. 2 1983, 79, 19. ( 8 ) MarEelja, S.;Radic, N. Chem. Phys. Lcrr. 1976,12, 129. (9) Gruen. D. W. R.; Marklja, S. J. Chem. Soc., Faraday Trans. 2 1983, 79, 225. Attard, P.; Batchelor, M. T. Chem. Phys. Lcrr. 1988, 149, 206. (IO) Kjellander, R. J . Chem. Soc. Faraday Trans. 2 1984, 80, 1323. (1 I ) Granfeldt, M. K.:Jansson, Bo.; Wennerstrbm, H. Mol. Phys. 1988, 64, 129.
its conclusions are irrelevant to the study of forces at small sep arations (where the repulsions are observed), these same authors later boldly suggested that the interaction of two confined dipolar layers across a solvent continuum is attractive over all separat i o n ~ . ' The ~ question of whether the image charge mechanism plays a primary role has yet to be settled and is still a current topic of debate.I4 Other molecular simulations have tested the water polarization model but have also failed to provide adequate support for that suggested m e c h a n i ~ m . ' ~ J ~ The collection of data on phospholipids suggests another molecular mechanism, related as much to properties of an entire amphiphile as to the whole bilayer. The mechanism proposed is, after the fact, an obvious one, made obvious by the conclusions of Rand and Parsegian' that the forces "do correlate in a systematic way with polar group identity and state of the hydrocarbon chain"! This statement, it tums out, is as much relevant to nonionic surfactant layers as to phospholipid layers. In short, there is definitely evidence to support an amphiphilic origin to the repulsions measured in these systems. One model that focuses on the properties of the entire bilayer was first proposed by Helfrich.'' Here the amphiphilic monolayers are treated as infinite undulating surfaces whose thermal fluctuations are sterically hindered in the presence of a second surface, thus causing a repulsive interaction. Unfortunately, the intrinsic separation dependence of this force (in the long wave limit) was a power law and not exponential. Israelachvili and Wennerstrom16have approached the problem on a molecular level. Their hypothesis is that the measured forces are essentially a result of steric interference between fingerlike projections of opposing bilayer surfaces. They derive a simple formula for the pressure, with an exponential-like behavior, having a dependence on the hydrocarbon/water interfacial free energy which, while being a continuum property, can give a consistent qualitative account of the trends in decay length and magnitude for different hydrocarbon chains, hydrated (neutral) headgroup, solvent, and temperature! The idea of a molecular origin directly involving the interaction between phospholipid molecules was indeed suggested earlier by Marra and Israelachvili." However, it was not considered a likely candidate2' as general opinion, based on an assortment of independent spectroscopic studies such as phosphorus 31 and deuterium NMR2*- and X-ray and neutron ~ c a t t e r i n g , *held ~ - ~ ~that the
'
(12) Attard, P.; Mitchell, D. J. Chem. Phys. Lcrr. 1987, 133, 347. (13) Attard, P.; Mitchell, D. J. J . Chem. Phys. 1988,88, 4391. (14) Leikin, S.;Kornyshev, A. A. J . Chem. Phys. 1990. 92,6890. (15) Gruen, D. W. R..; MarEelja, S.; Pailthorpc; B. A. Chem. Phys. Lcrr. 1981, 82, 315. (16) Kjellander, R.; MarEelja, S. Chem. Scr. 1985, 25, 73. (1 7) Helfrich, W. Z . Naturforsch. 1978, H A , 305. (18) Israelachvili, J. N.; Wcnnerstrbm, H. Lungmuir 1990, 6, 873. (19) At the time of writing this article we were made aware that simula-
tions of a system similar to ours had been carried out and preliminary mdts submitted elsewhere (Nilsson, U.; Jansaon, Be.; Wenncrstram, H. Faraday Discuss. Chem. SM., in press). Their model has the same physical basis and their findings are also similar. (20) Marra, J.; Israelachvili, J. N . Biochemistry 1985, 24,4608. Scc also: Marra, J. J . Colloid Interface Sei. 1985, 107, 446; 1986, 109, 11. (21) Gruen, D. W. R.; Marklja, S.; Parsegian, V. A. Water Structure Near The Membrane Surface. In Cell Surface Phenomena; Marcel-Dekker: New York, 1984; Section IV, pp 51-91. (22) Seelig, J.; Gally, H . 4 . Biochemistry 1976, 15, 5199. (23) Seelig, J.; Gally, H.-U.; Wohlgemuth, R. Biochim. Biophys. Acra
1977,467, 109. (24) Seelig, J. Biochim. Biophys. Acra 1978, 515, 105. (25) Seelig, J.; Seelig, A. Q.Rev. Biophys. 1980, 13, 19. (26) Ycagle, P. L.; Hutton, W. C.; Huang, C.; Martin, R. 9. Biochemistry 1977, 16,4344. (27) Yeagle, P. L. Acc. Chem. Res. 1978, 1 1 , 321. (28) Parson. R. H.: Pascher. I. Nature 1979. 281. 499. i29j Phi!lips,' M. C.; Finer, E. 0.;Hauser, H. Bi6chim. Biophys. Acta 1972, 290, ~ 9 1 . (30) Levine, Y. K.; Bailey, A. 1.; Wilkins, M. H. F. N o r m 1968,220,577. (31) Bllldt. G.; Gally, H.-U.; Seclig, J.; Seclig, A. Nature 1978,271, 182. (32) Bllldt, G.; Gally, H . 4 . ; Seelig, J.; Zaccai, G. J . Mol. Biol. 1979,131, 673.
Simulation of Flexible Zwitterions dipolar head of phospholipids (directed from the phosphate to the nitrogen) lies 'almost parallel to the plane of the bilayer"! This orientation was argued to persist over a wide range of conditions. The Israelachvili and Wennerstr6m hypothesis is now open to confirmation by simulation and is the subject of this communication.19 We aim to show that repulsive forces between the flexible zwitterionic layers are indeed a natural consequence of the partial diffusive freedom of the molecules and that these forces have a magnitude and separation dependence very similar to those measured. We include the flexibility of the amphiphilic molecule in these simulations by extending a previous model" in two respects; neither the distance between the ionic centers nor the location of the zwitterion with respect to the interface is fixed. Although simple, the model shows a definite and consistent connection between the state of the 'hydrocarbon core", the 'hydrated" size of the charge groups, and the magnitude of the force and its decay constant, as found in experiments and in much the same manner as proposed by Israelachvili and Wennerstram.I* At the same time, by comparing our simulated 'headgroup" conformations with experiments, we show that the interpreted near-parallel alignment is consistent with the model and by no means denies the possibility of an entropic repulsion.
Description of the Model (a) The Physical System. Although now over a decade old, Tan ford's monograph34 still provides a concise and informative picture of general amphiphilic behavior. A review of some principal features of lipids is useful for the development of a suitable model molecule. The fluidity of lipid bilayers is strongly correlated to length and type of hydrocarbon chain present. For saturated chains, the shorter the chain length, the lower the transition temperature, T,, from an ordered to a disordered bilayer state. For unsaturated chains the transition temperature is substantially less than for saturated chains of the same length and less still for higher degrees of unsaturation. Such correlation is readily evidenced by the diffusion rates of lipids in the plane of the bilayer which can change by 2-3 orders of magnitude as one crosses Tc.35 Even so, at room temperature the diffusion constant in the case of saturated phospholipid bilayers ( T , = 42 "C) is quite high, of the order of IO-* cm2/s. Though the interior may be fluid, the time scale for the movement or exchange of lipids either out of the bilayer into solution or between the two sides of the bilayer has been quoted in terms of hours or even days. The first of the latter processes calls for the transfer of an entire hydrocarbon chain into aqueous solution which, in the case of naturally occurring phospholipid systems containing generally low free lipid concentrations (a cmc of approximately M), has associated with it an unfavorable free energy of transfer of about 15 kcal/mol. The second process also requires an unfavorable event, the passage of the hydrophilic headgroup through the hydrophobic core. For phospholipids, the requisite activation energy may be estimated by the electrostatic free energy of transfer as provided by the Born model of ion solvation.2J6 A value of about 40 kcal/mol per charge can be expected for ions with an approximate radius of 2 Furthermore, for zwitterions there is an additional transfer cost of the order of e2( I/tk - l/tw)/R, where R is the (assumed) fixed distance between the two contributing ions. A fluidlike membrane, however, has a further consequence which has indeed been explicitly considered Though (33) BOldt, G.; Seelig, J. Biochemistry 1980, 19, 6170. (34) Tanford, C. The Hydrophobic Effect: Formation of Micelles and Biological Membranes, 2nd ed.; Wiley-lnterscience: New York, 1980. (35) Yeagle, P. The Membranes of Cells; Academic Press: New York, 1987. (36) Bockris, J. OM.; Reddy, A. K. N . Modern Electrochemistry; Plenum Press: New York, 1970 Vol. I . (37) Rashin. A. A.; Honig, B. J . Phys. Chem. 1985, 89, 5588. (38) Aniansson, G.E. A. J . Phys. Chem. 1978,82, 2805. (39) Gruen, D.W. R.; de Lacey, E. H. B. In Surfactants in Solufion; Mittal, K.L.,Lindman, B., Eds.; Plenum Press: New York, 1984: Vol I , pp 279-306.
The Journal of Physical Chemistry, Vol. 95, No. 16, 1991 6353 h
2
-21
2
4
6
8
Distance (A) Figure 1. A plot of the pair probability distribution (in arbitrary units) of carboxylic oxygens in glutamic acid obtained from a molecular dynamics simulation of a small peptide in aqueous solution. In this plot the value of rOsis given as the position of the peak in the probability distribution.
the exchange of phospholipids between the bilayer and solution is unfavorable a t physiological concentrations, limited spatial fluctuations of the amphiphilic ion, with a component in a direction normal to the average hydrocarbon/solvent interface, is permissible. Such random movement in a third dimension is entropically favorable a t a relatively small cost of about 400-800 cal/mol per CH2 group exposed to an aqueous envir~nment,'~Le., in unit distances of approximately 2 This anticipated spatial fluctuation is likely to be fast and, together with some residual rotational ability along the hydrophilic chain, shows up in X-ray and neutron scattering studies as a smooth d i s t r i b u t i ~ nleaving ~,~~ a surface which Tanford calls locally 'rough ...containing fluctuating peaks and troughs". Naturally, this degree of spatial fluctuation is strongly dependent on the nature of the hydrocarbon chain, the temperature, and the particular solvent environment. As to the zwitterionic headgroup itself, the linear chain linking the cationic group to the anionic group, while quite short (in phosphatidylcholine (PC) and ethanolamine (PE), two methyl groups and an oxygen between the phosphorus and nitrogen groups), retains some rotational ability resulting in an internal conformational distribution. The most probable combination of torsion angles is determined from a balance of electrostatic attraction, steric repulsion, and conformation energy considerations. One expects a tight but non-delta-like distribution of one charge from the other, much the same as is shown in Figure 1 (obtained from a molecular dynamics simulation of a small peptide46). (b) The Theoretical Model. Given that the hydrophobic moiety of lipids is understood to be the primary mechanism driving agg r e g a t i ~ nit, ~is~important that the intimate connection between the headgroup to which we focus attention here and the hydrophobic tail not be overlooked, as any aggregate system is only partially dependent on properties of either part of the molecule. While it is difficult to include the entire character of the hydrophobic chain, its role must be considered in some way. On the basis of the above review we can formulate a simple model of a zwitterionic (mono)layer. Firstly, each zwitterionic ~
(40) Gruen, D. W. R. J . Phys. Chem. 1985,89, 146. (41) Gruen, D. W. R. J . Phys. Chem. 1985,89, 153. (42) Karaboni, S.; OConnell, J. P. J . Phys. Chem. 1990, 91, 2624. (43) Karaboni, S.; O'Connell, J. P. Lungmuir 1990, 6, 905. (44) Stigter, D.; Dill, K. A. Lungmuir 1986, 2, 791. (45) These values are very much dependent on the nature of the headgroup, chain length and type, solvent, and temperature. The values also vary, depending on whether the lipids are removed from a micellar aggregate (in which case the headgroups are always in polar solution) or from bulk hydrocarbon. For a particular set of conditions, the quoted energy cost per CHI group is based on an asymptotic linear least-squares fitting of the free energy transfer curve plotted as a function of chain length. However, it is recognized that the exposure of the methylene groups closest to the headgroup is less costly. Thus, for short exposure distances the functional form is expected to differ from linear dependence. (46) Ulner, M. Private communication.
Granfeldt and Miklavic
6354 The Journal of Physical Chemistry, Vol. 95, No. 16, 1991
head, with a chemical composition of the form A+(CH2),B- (e.g., PC has A+ = N+(CH& n = 2, and B- = PO3-), is modeled as two charges, *e, embedded in hard spheres of radii a,. These two spheres are then connected by means of a simple harmonic potential uh(r) = Kh(r - rq)2
(1)
acting between their centers, a distance of r apart. Kh (subscript denotes headgroup) is some effective "spring" constant whose magnitude is a measure of the distribution of one charge from the other, about an equilibrium center of oscillation, rq (see Figure 1). These constants may be readily set for a particular system either by me; 11s of comparison with independent experimental data4' or via a more detailed simulation of the specific molecule (as used to generate Figure 1). Both Kh and rq possess some intrinsic knowledge of the internal structure of the molecule. In application to phospholipids the headgroup chain is marginally flexible with the equilibrium separation, rq, not too dissimilar from the fully extended, all-trans, length. We can therefore set the parameters Kh and rq to give a realistic average charge-charge distance in the range 4.5-5.5 A, e.g., for a comparison with the phospholipids PC and PE. Two of the examples studied in this paper, Kh = 0.0546 J/m2, rq = 0.0 A; Kh = 0.437 J/m2, res = 5.0 A, satisfy this criterion. Assigning appropriate values to the sizes of the positive and negative particles, a*, is made difficult by the unknown extent to which each ion is hydrated; such contributions are not easy to transpose from one experimental situation to another. In this paper we shall first set the negative/positive charge groups to have equal 2.0-A radii and then examine the effect of increasing the size of the positive ion. A second potential is needed to describe the way in which the "headgroups" are connected to one or more hydrocarbon chains. Our discussion above suggests a model potential acting at the anionic center which is an increasing function of distance from the average hydrocarbon/solvent interface. We have chosen, primarily for numerical convenience, another harmonic potential u,(r) = Kg2
(2)
between the anion and the grafting point, although it should be immediately clear that beyond a distance equal to Imax, the maximum length of the hydrocarbon chain, for which the amphiphile is presumably completely removed from the bilayer, the potential should in principle be constant. Experimental studies have shown that the cost of removing an amphiphile either from pure alkane or from a micelle is asymptotically a linearly increasing function of n, the number of methylene groups.34 This fact has been utilized to estimate the energy cost of perpendicular motion of an amphiphile out of an aggregate.41*42*44*45 An estimate for the parameter Kc (subscript denotes chain) can be calculated simply by equating the measured total free energy cost of the complete removal of a phospholipid from an aggregate into solution (-15.1 k ~ a l / m o l ~to~ )the product Kclmax2.. For example, suppose that for the phospholipid dipalmitoylphosphatidylcholine (n = 15) I,,, is given by the all-trans chain length of 18.3 A; then K, would have a value of 0.032 J/m2. Even allowing for some shrinkage of the chain by some 20% (from its self-association to reduce the area of exposed hydrocarbon), we get a value of K, = 0.049 J/m2; with 40% shrinkage, K, = 0.087 J/m2. It is unfortunate that the extremely low cmc of naturally occurring phospholipids all but precludes evaluations of the free energy of transfer across a range of chain lengths and/or headgroups. However, the above values for K, give reasonable estimates and are therefore a useful guide for a parameter study. Three further features of the model system should be made explicit. First, the aqueous environment in which the ions are assumed immersed is treated at the continuum level. That is, a dielectric of constant relative permittivity, ew, is used as substitute, and no explicit involvement of the solvent is considered. Ion solvation effects therefore only show up as effective ionic sizes, (47) E.g.: Chevalier,
Y .; Perchec, P. L. J. Phys. Chem. 1990, 94, 1768.
Figure 2. A model zwitterion molecule with connecting "springs". The particle labeling of "0" for connecting point, "1" for negative group, and "2" for positive group is shown. Interparticle distances calculated in the simulation are also defined.
Secondly, as phospholipid membranes do not readily solvate ions, our effective "hydrocarbon/water" interfaces at z = k b will be treated as impenetrable hard walls. Finally, the contribution from image charges, which was studied in our previous work, has not been included. One reason is the increased computational effort involved; another is that in those investigations, with zwitterions located on the water side of the interface, we have already shown that at 2.0 A into the water phase the image charge effect is weak.' A schematic representation of our modeled headgroup is shown in Figure 2 along with a suggested particle labeling referred to in the text.
a,.
Simulation Details and Calculated Properties Monte Carlo Simulation. The extremely low cmc value of naturally occurring phospholipids allows one to study the total system of two interacting layers of zwitterions, with a given number of lipids per unit area, regarded as virtually isolated from a bulk lipid solution. Our calculations are based on the canonical Monte Carlo technique using a Metropolis algorithm.48 The Monte Carlo box was a rectangular prism, and periodic boundary conditions were applied in the x,y directions. During the course of the simulation all potentials acting between the hard sphere and headgroup particles confined to the simulation box are evaluated. The total energy within a single zwitterion, Uintra, includes the electrostatic and both harmonic potentials Uintra
= Kc(Ir1 - r01)2 + Kh(lr2-
rll
- rq)2 - e2/4w$w(r2
- rll (3)
where r1,2are the vector positions of the charge spheres (see Figure 2). ro locates the point of connection of the zwitterion to its surface. Only electrostatic potentials act between the zwitterions so that, provided there is no hard-sphere overlap, the total intermolecular energy is izj
-
(4)
All particles are given random displacements, including the points of connection which only move in the plane of the surface. Apart from explicit interactions evaluated inside, all charged particles experience an external potential arising from the zwitterions present outside in the infinite planar system. The action of the harmonic potentials results in distinct, nonuniform distributions for each ion type (see Figure 3). These distributions combine to generate a nonzero average electrostatic field which acts on the ions within the box. The external potential, calculated within the mean-field approximation, is updated with each new simulation run, taking the current concentration profile found within the box as the source of the externally applied field in the (48) Metropolis, N.A,; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A.; Teller, E. J. Chem. Phys. 1953, 21, 1087.
The Journal of Physical Chemistry, Vol. 95, No. 16, 1991 6355
Simulation of Flexible Zwitterions subsequent run. Self-consistency is normally achieved within two
or three of such iterations steps. The simulations were performed with a total box content of 264 particles: 88 anions, 88 cations, and 88 connecting points, giving 44 zwitterions per surface dispersed at a density of 1 in every 60 A2. Equilibration runs of 75 000 moves per particle or
225 000 moves per zwitterion were performed prior to all production runs which in turn were of the order of 150-200000 moves per particle. Individual displacement parameters were assigned to the different particles (0, 1, and 2) and varied to give a reliable acceptance rate. The temperature of the system was kept constant at T = 298 K, and the dielectric constant was that of water, cw = 78.3. Conformational Properties. Properties of the system are calculated at regular intervals in the simulation and averaged upon completion. Internal zwitterion statistics are accumulated on the average positive-negative separation, r I 2 (averaged over all molecules), the distance of each particle from the wall to which they are connected, zIoand z2@and finally the average difference in z value of the two charge particles, zI2.These are defined generally by the formula
where [ denotes the configuration property, r or z, s, t = 0, 1, or 2, and the sum is over all molecules in the box, Nmil,.In the case of zI2no absolute value is taken. For convenience, these variables are descriptively defined in Figure 2. Each average is accompanied by a variance, or m n d moment of the average, u, which measures the width of the distribution about the mean, evaluated in the usual way:
Besides the internal averages, the density distributions of particles 1 and 2, pl(z) and &), respectively, were calculated. In general, the calculated averages and distribution functions were quite reliable, exhibiting little or no variation between iteration runs. Osmotic Pressure. The most important quantity we study in this system is the force per unit area between the two opposing zwitterionic layers. Within the primitive model, i.e., with a continuum representation for the solvent and treating the particles as charged hard spheres interacting via pair potentials, the internal osmotic pressure can be evaluated by means of a contact equation
A similar expression was derived and used in the context of polyelectrolyte solutions between charged planar walls49*50 and appars as a straightforward generalization of the well-known contact condition for simple electrolyte system^.^^*^^ In eq 7, pl(2)(bals) is the density at closest approach to the wall, z = b, of particle 1 (2) and u, is the grafting potential defined earlier; &,/ab is the component of the corresponding stretch force in the direction normal to the walls, and the summation concerns only those molecules (of number NZwil1/2)connected to the wall at z = b. For flexible zwitterionic particles connected by “spring” potentials this pressure relation takes the form49 1 “ttll2 Pint = Pi(b - a J k T - 2K,- E (6 - z i ) (8) 1-12
A
j-1
(49) viklavic, S.J.; Woodward, C. E. J . Chem. Phys. 1990, 93, 1369. (50) A k a y n , T.;Woodward, C. E.; Jansson, B.J . C h p . Phys. 1989,91, 2461. Miklavic, S. J.; Woodward, C. E.; Jbnsson, B.;Akesson, T. Macromolecules 1990, 23, 4149. (51) Henderson, D.;Blum, L.; Lebowitz, J. L. J . Electround. Chem. Interfacial Electrochem. 1979, 102, 31 5 . (52) Janmn, B.; Linee, P.; Akesmn, T.; Wennerstram. H. In SurJucrcmrs in Solution; Mittal, K. L., Lindman, B.. Eds.; Plenum: New York, I984 Vol. 3.
Although such a simple form for the pressure relation is available, it is still numerically one of the most difficult thermodynamic properties to evaluate. In some cases a great deal of accuracy is required to determine the contact density since the difference in the two terms in eq 8 can be smaller than the statistical error. Naturally, the greatest uncertainty occurs under conditions of near-zero pressure, that is, around crossover points from attraction to repulsion and generally at large separations (normally 2 4 0 A). In principle, the accuracy can be increased by simply increasing the length of the simulations. The second term on the right-hand side can be evaluated with quite high accuracy under most conditions, but the term involving the particles’ contact density presents a limitation to the accuracy around these regions of net zero pressure. Fortunately, in the greater part of separations of our current interest where experiment has found extremely large repulsive pressures (between 5 and 30 A), we find that the contact density and therefore the pressure are much larger than the statistical error.
Results and Discussion The behavior of this model system of two opposing zwitterionic layers can best be understood by examining parameter space and in limiting cases making associations with well-known systems. We focus attention on the obvious “unknowns”, the parameters in the connecting potentials, discussing two cases; variable K, = Kh = K with rCs= 0.0, and fixed Kh = 0.437 J/m2, rCs= 5.0A, and variable K,. In the former case variations have also been made in headgroup size and charge. We will present the results in two sections, beginning with (a) conformational properties followed by (b) interlayer interaction. (a) ConformationalProperties. (i) Equal Spring Constant,K. Simple electrolyte and neutral hard-sphere solution^^^"^ are convenient states of reference for our discussion, the reason being that electrostatic and hard-core interactions govern this system’s behavior in a similar way. Let us begin with the simplest case: the two charges on each molecule are connected by the potential of eq 1, with rOs= 0.0, and the two spring constants set equal. We can study the competing effects of electrostatics and hard-core interactions as the parameter K is varied. In setting rq equal to zero, the particles have a minimum separation limit set by the sum of their hard-sphere radii, in this case 4 A. Figure 3 shows the density distributions of the two ions that result from the variation of K: the distribution of particle 1 is shown in Figure 3a and that of particle 2 in Figure 3b. For a sufficiently small value of K = 0.0002 J/m2, all particles act independently, leaving undefined the internal molecular structure and any association with a given wall. The result is a hard-sphere monovalent electrolyte which acts indiscriminately of ion charge and therefore shows identical profiles for the positive and negative particle distributions. The fact that the profiles are not uniform is a consequence of hard-sphere effects. At this separation, which is 3.75 times the hard-sphere diameter, there is rmm enough for three distinct layers of the particles: two intense layers at the planes of closest approach and a more diffuse layer in the center. This is the preferred packing arrangement at this volume fraction (-0.1).53 As K is increased, the internal molecular structure starts to be defined and an organizational process begins. For a K value of 0.01 18 J/m2 negative particles are pulled in to their associated wall (Figure 3a), forcing some of the positive particles outpredominantly a volume exclusion effect (Figure 3b); the ejected positive ions, still relatively free, then seek their own preferred packing arrangement midway between the walls increasing the density there. At still higher values of K = 0.0546 J/m2, both particles are pulled in, e.g., the density of anions at the wall increases 2-fold, achieving a region of intense negative charge, while diminishing substantially at the center. The tight springs increase the likelihood of steric repulsion between the positive (53) Sarman, S.;Kjellandcr, R. Chem. Phys. Lerr. 1988, 149, 102. (54)Kjellandcr, R.;MarEelja, S. J . Chem. Phys. 1988, 88, 7138. (55) Svensson, B.;Jansson, B.;Woodward, C. E. J . Phys. Chem. 1990,94, 2105.
6356 The Journal of Physical Chemistry, Vol. 95, No. 16, 1991
Granfeldt and Miklavic
C
2
4
-5
-4 -3 -2 -1
2
1
0
3
5
4
Distance (A)
-5
-4 -3 -2 -1 0 1 2 Distance (A)
3
4
5
Figure 4. Density distribution of negative (1) and positive (2) particles for the case K, = Kh = 0.0546 J/m2 and rOs= 0.0. The curves of small (large) magnitude at the center correspond to 1 (2) particles in the charged and uncharged states. Full lines refer to particles in the charged state, and dashed lines refer to particles in the uncharged state. The surface separation is 15 A.
TABLE I: Internal Structural Distances, Averaged over the Total Number of Zwitterion Molecules, for the Parameter Set K, = K, = 0.0546 J/m2. r.. = 0.0" av uncharged charged r(1-2), A 4.8 (0.8) 4.7 (0.7) 3.3 (1.1) z(o-I), A 3.1 (0.9) 5.1 (1.9) 4.8 (2.0) z(0-21, A z(l-2), A 2.0 (1.9) 1.5 (2.0) 18.4 9, deg 24.7 ~~
~
& & 2
-5
-4 -3 -2 - 1
0
1
2
3
4
5
Distance (A) Figure 3. Density distribution of particles for the spring parameter set Kc Kb = 0.0. Part a shows the negative particle density for the casLydL%0002 J/m (dot-dashed line), K = 0.0118 J/m2 (dashed line), and K = 0.0546 J/m2 (solid line). Part b shows the positive particle distribution under corresponding conditions. The surface separation is 15 A. particles and the negative particles present near the wall, striving to expel the former. This directly competes with the electrostatic attraction between the positive charges and the built-up anionic layer, drawing them in toward the surface, and with the springs between a pair of particles keeping them together. The total combination results in a splitting of the central peak in the positive charge distribution! This is predominantly due to the action of the spring itself while electrostatic interactions modify the m a g nitude of this effect. It is clear that each of the two new peaks is associated with a given surface. At this surface separation the peaks are separated by a distance of the order of the radius of the spheres and not their diameter, which may suggest that these ions on average sit interstitially in the middle, with the ion centers of one layer in the plane of the perimeter of the other layer. In Figure 4, the densities for the two particles (1 and 2), at K = 0.0546 J/m2, are shown both with and without charges. From the figure it is clear that the concentration of the negative ions close to the walls is depleted as a result of electrostatic repulsion within each anionic layer, while the attractive interaction between this layer and the positive charges associated with the same surface draws the latter in, increasing the density by a factor of 1.5 from the uncharged case. Here it is quite evident that electrostatics only modifies the already split central peaks in the particle 2 distribution. The changes in the distribution functions are also reflected in changes in the average internal distances given in Table I. The response to the introduction of charge is as one would expect. The depletion of negative charges near the wall is reflected by a slight increase in the average distance of particle 1 from the wall, z , ~ , Similarly, the movement of particle 2 closer to the wall is reflected in decreases in r I 2and zO2. The value of zO2quoted in Table I is effectively an average over the two peaks per surface; thus, the difference reflects the variation in relative population at these peaks rather than their change in position. It-isclear from the figure and the table that electrostatic interactions reduce the average
'The list demonstrates the dependence on the state of charge at a separation of 25 A. The labels 0, I , and 2 correspond to the labeling in Figure 2. The values in parentheses are standard deviations, i.e., variances of averages as calculated with eq 6. orientation of the zwitterion with respect to the plane of the bilayer, ), is given measured crudely by the angle qj = sin-' ( z I 2 / r l 2which in the last row of Table I. A uniform orientation distribution gives an angle of 30° to the plane. We found this in the case of charged zwitterions fixed a certain distance from the interface.'' (This is the same as for the uncharged case here.) For charged molecules, we now find them more likely to have a relative orientation more parallel to the plane, for reasons discussed above. In order to obtain a smaller angle with electrostatic interactions, it is thus necessary for either the zwitterion to be allowed out-of-plane fluctuations escaping intralayer steric constraints, in which case the K, value of 0.0546 J/mZ evidently allows for sufficient movement, or the hard-core radii to be reduced.36 Using the same value of 0.0546 J/m2 for the two spring constants, we have also investigated the effect of headgroup size. This is of interest as the difference between PE and PC in the present model is a larger hard-core radius of the positive ion for PC. We compare the case of particles 1 and 2 with radii of 2.0 A with that where the radius of particle 2 has been increased to 3.0 A. Figure 5a shows the negative particle density which is hardly affected. For the positive particle (Figure 5b), however, the change is quite large. The plane of closest approach increases by 1 A, reducing the available volume and thus causing an overall increase in the density. The off-surface peak is also moved by approximately the same amount. (ii) Variable K , for Fixed K h and Finite r In this case we vary the spring constant K, at fixed Kh = 0.4% J/m2 and with an equilibrium oscillation distance, rq, of 5.0 A. Table I1 summarizes the average distances r I 2 ,zl0, zlz, and zo2;r I 2of course remains constant with K, while z12diminishes indicating that, in the limit of zero K,, the vector directed between the connected particles is equally likely to be pointing up as down. The average distances of the particles from their connecting wall, zol and zoz, approach a value close to half the surface separation, as expected.
.
(56) Jonsson, B.; Attard, P.; Mitchell, D. J. J . Phys. Chem. 1988, 92, 5001.
Simulation of Flexible Zwitterions
The Journal of Physical Chemistry, Vol. 95, No. 16, 1991 6357
-55
1.0
.e,
t,
a
.-P ;0.5 3 2
g
.-x 0.0 -8
-6
-4
-2
0 Distance
2
4
6
8
(A)
h
3
!i
p - 8 - 6 - 4 - 2
0
Distance
2
4
8
6
(A)
Figure 5. Density distribution of particles for the spring parameter set KE= Kh = 0.0546 J/mz. Full lines refer to the case 2.0/2.0-A radii of negative/paitivc particles; the case of 2.0/3.0 A is represented by dashed lines. Part a (b) shows the negative (positive) particle density. The surface separation is 20 A.
TABLE II: Dcpendme of Structunl Averages at Fixed K, = 0.437 I,, = 5.0 A on the Spring Constant K,O Kc = Kc = Kc = Kc =
J/m2 a d
0.437
av
J/m2
r(I-2),
A z(&I), A ~(&2), A 2(1-2), A
5.1 (0.7)
4, deg
2.2 (0.2) 4.5 (1.6) 2.3 (1.6) 26.3
0.055 J/m2
0.023
5.1 (0.6) 3.3 (1.1)
5.1 (0.6)
4.9 (2.1) 1.7 (2.1) 18.9
0.0002
J/m2 4.0 (1.6) 5.2 (2.4) 1.2 (2.3) 13.5
J/m2 5.1 7.5 8.2 0.8 8.8
(0.6) (3.4) (-) (2.6)
OThe separation of the surfaces is 15 A, and the labels 0, I , and 2 correspond to the labeling in Figure 2. Again we see the significant effect of the fluctuations on the angle 4. For a relatively strong spring of K, = 0.437 J/m2 the average is again similar to the value of 30" found with fixed zwitterions.'l Distribution functions for this case are not shown as they appear similar to case i except for obvious modifications related to the slightly extended distance between negative and positive particles. Finally, we have only presented data a t certain separations; calculations performed a t others show that the separation dependence of these internal averages and orientations is very small. (iii) Relation to Experiments. The spectroscopic methods deuterium (D) and phosphorus31 (31P)NMR,22-25J7J8 NMR utilizing the nuclear Overhauser effect,265' X-ray crystallography,28 and neutron ~cattering~l-'~ have provided most of the data on conformational structure of phospholipids in bilayer membranes, particularly PC and PE. The general conclusion reached in these studies specifies that the PC and PE headgroup conformation, given in terms of the vector directed from the phosphor atom to the nitrogen (denoted by PN),lies "parallel" or "almost parallel" (57) (58)
Skarjunc, R.; Oldfield, E. Biochemistry 1979, 18, 5903. Skarjune, R.; Oldfield, E. Biochemistry 1982, 21, 3154.
to the average plane of the bilayer. To test our results against this experimentally defined condition, it is first necessary to quantify on each experimental account, if possible, the meaning and significance of "almost parallel". The X-ray structure of DMPC dihydrate crystals (kept between 10 and 15 oC)28shows the headgroups to be alternating between two average orientations: the PN vectors of length 4.3 and 4.5 A making angles of 17 and 27" with the plane of the bilayer. These values (termed "almost parallel" by Seelig and SeeligZ5) can be compared with the values listed in column 2 of Table I and Table 11, bearing in mind that the simulations were performed at a higher temperature (25 "C) and our quoted angle (like the values of zI2)are averaged over the two peaks per surface in the positive particle distribution (Figure 4). As stated in the description of the theoretical model, we feel that these results represent a fair comparison with the real system. From D and 31P N M R studies of selectively deuterated PE22 and PC23*24 a set of intramolecular torsion angles were deduced by theoretically reproducing the experimentally found maximum chemical shift anisotropy. This set of angles specifies a headgroup orientation parallel to the membrane. Skarjune and Oldfield,s7 however, have given criticism to this conclusion, suggesting that orientation space had not been properly sampled and that, in fact, substantial regions of this phase space were found to be consistent with the observed experimental data while not always giving parallel alignment. This has placed some doubt on the quantitative information previously presented, and the issue, all but ignored,2s has not been resolved to the authors' k n ~ w l e d g e . ~ ~ * ~ ~ By exploiting the nuclear Overhauser effect in N M R and following a logical experimental line, Yeagle et reason that the PC dipole is predominantly oriented "approximately parallel" to the bilayer as a result of intermolecular electrostatic interactions between the positive N-methyl group and negative phosphates of neighboring molecules in the same monolayer. Furthermore, their findings indicate that with the addition of cholesterol, as a spacer molecule to increase the average separation between the phospholipids and thus reduce intermolecular interactions, the headgroup orientational freedom is increased. (This incidentally is contrary to earlier claims of orientational indifference to the addition of c h o l e s t e r 0 1 . ~ ~ * ~While - ~ ~ ) still only qualitative, their findings do indicate that the PN vector is more likely directed at small angles to the bilayer rather than perpendicular, depending on the strength of intermolecular electrostatic interactions. Although we have not here presented any data appropriate to the incorporation of neutral spacer molecules, the results given in Figure 4 and Table I, where the effect of charge is demonstrated, do give some confirmation to the conclusions of Yeagle et al. More quantitative information on the preferred orientation has come from the series of neutron scattering studies of selectively deuterated In these, average distances (and their variances) from the center of the hydrocarbon core of selectively deuterated groups have been determined by theoretical modeling of corrected scattering data. As more information is available on the PC lipid,3'*32the discussion and comparison to follow are based on this data. The experimental results indeed indicate that, while significant changes in the state of the hydrocarbon chain do occur with either increased temperature or degree of hydration, the average orientation of the headgroup does not show as great a response: the difference between the distance from the center of the hydrocarbon core to the nitrogen group and to the first methylene group above the phosphate varies from 0.6 to 1.5 A. Given that the phosphate itself (with the oxygens sharing the negative charge) lies somewhere between this methylene and the glyceryl entity (a region of estimated width between 1.4 and 3.6 A, depending on conditions), it is not unreasonable to expect an average difference in distance between the charges ranging from approximately 1.3 to 2.4 A, with increased hydration and temperature. These values have a comfortable correspondence (59) Dill, K. A.; Stigter, D.Biochemistry 1988, 27, 3446. (60) Franks, N. P. J . Mol. Biol. 1976, 100, 345. (61) Worcester, D. L.; Franks, N. P. J . Mol. B o / . 1976, 100, 359.
Granfeldt and Miklavic
6358 The Journal of Physical Chemistry, Vol. 95, No. 16, 199'1
with our simulated z12values given in the second column of Tables I and 11. Biildt et al.32have also determined the variance about the mean in one case; for highly oriented samples a t low water content, their calculated spread (assuming a Gaussian distribution) of one deuterated methylene group halfway along the head is given a value of 2.12 A at 20 OC, increasing to 2.55 A at 70 "C. Unfortunately, we can only quote changes with separation: our two representative systems have C T , ; variances ~~ increasing up to 2 A, with separation, certainly below the range determined by Biildt et al.32 Taking into account the variance about the mean position of the topmost carbon on one hydrocarbon chain (0.9 A), Biildt et al.32concluded that the above spread in the distribution of the headgroup methylene is partly due to the motion of the whole molecule, in the direction perpendicular to the bilayer plane. On the basis of our own calculations we can support their conclusion. To this we may add our own interpretation that the average orientation of the PN vector, lying close to the bilayer plane, is itselfa consequence of this motion, with electrostatic interactions adding an extra contribution. Stigter and Di11,59u62however, have fully accepted the (near) parallel orientation as a consequence of a totally different phenomenon. They have proposed a model in which the tail of the dipole (the negative phosphate group) sits at the interface and its head (the positive nitrogen group) is directed into the hydrocarbon phase of the membrane (and going further into the hydrocarbon with increasing temperature) due to the hydrophobic affinity of the methyl cage surrounding the (charged) introgen for the hydrocarbon region. They too have provided positional averages and variances which, by and large, fall within experimental error. While we cannot offer any direct evidence of our own to refute their model, the neutron scattering results31J3suggest one reason why it is not feasible: the headgroup deuterated carbons are known to be about 5-6 A a h e the average plane of the hydrocarbon/water interface with a glycerol/water region intervening. Furthermore, this distance increases along the length of the hydrophilic moiety. In conclusion, the results of our model by no means contradict the experimental findings concerning the headgroup orientation. However, it obviously does not lend itself to support the interpretation of a rigidly parallel alignment of the headgroup. (b) Interlayer Interaction. As the numerical value of the total pressure is the same at any plane between the two walls, it is often informative to discuss its parameter dependence in terms of the different contributions at two particular planes. At either wall, the pressure is given by eq 8. At the center, the pressure expression has a far more complicated form with contributions from entropy (p(O)kT), hard-core interaction, and electrostatic and harmonic correlations across the midplane. (For a discussion on the latter see, e.g., ref 50.) As the dominant repulsive contribution to the interaction is measured by the concentration at the midplane, it is obvious that a distribution of particles arising from the movement of the molecules normal to the plane of the layer is an important aspect to consider. Once there is an overlap in the particle distributions from the opposing layers a repulsion is possible. In Figure 6 we see the pressure for varying K , = Kh = K or K , at fixed Kh for a surface separation of 15 A. The curves tend at zero spring constant to the limits of a simple electrolyte solution and dipolar solution, respectively, for the two cases mentioned above. Connecting the charge particles into "free" zwitterions (the zero K , limit with fixed Kh) reduces the pressure compared with the simple free electrolyte case. This is caused by a lower contact density as a result of regions of configuration space being ~
~~~~~
(62) Stigter, D.; Dill, K. A. Lungmuir 1988, 4, 200. (63) We attempted to employ the model of ref 1 1 , excluding image charge
effects, at smaller surface separations than studied previously. The pressure could be calculated by using a contact equation (see ref 49), only if the size of the negative ion was set to zero placing it on the plane of the hard wall. However, this state makes for a difficult comparison with the model introducsd here. As we are missing one repulsive contribution, the total repulsive pressure plotted as open circles in Figure 8 would undoubtedly be extended to much larger separations and be of greater magnitude had the negative ion been of finite size.
0
-4
-3
-2
-1
0
Log K, Figure 6. Dependence of the pressure on the spring constant, K,. The filled squares represent the case of equal, variable springs, K, = Kh and Tal = 0.0 A. The open squares represent the case of fixed rq = 5.0 A and Kb = 0.437 J/m2 and variable K,. The surfaces are separated by 15
A; the particles have equal hard-core radii of 2 A. The pressure is plotted on a density scale in molar units. The conversion factor to dyn/cm2 is 2.478 X 10'.
excluded to the "orientable" molecules by the wall. As the spring constant(s) is (are) increased, the two curves converge as now the charge pair in the case of variable K also become coupled. For K , greater than 0.0015 J/m2, the springs attached to the walls draw the particles in to the surfaces and consequently the attractive stretch term of eq 8 begins to contribute to the pressure. At the midplane, the result is a reduced density (and also reduced hard-core interaction); the contribution from the charge correlations across the midplane is then relatively greater. The curve for the case of fixed Kh and rq = 5.0 levels off a t a slightly lower value, again presumably as the ions within a molecule are separated by a larger distance and therefore have a greater configurational constriction near the wall. For all these values of spring constant a repulsion results. How then does the force depend on separation? Characteristic of all the measured forces in phospholipid systems is an exponential repulsion with a decay length ca. 2 A, sometimes affiliated with the size of the water molecule. A second characteristic of these systems is that the bilayers swell up to a finite distance of d,,,, in excess water. Most of the postulated mecha n i s m ~ ~ ,have ~ * ~an ~ -exponential ~* separation dependence of some description; by adding an attractive van der Waals potential one can also account for the finite lamella swelling distance, d-. A closer inspection of the composite experimental data reveals a far more intricate situation involving the hydrocarbon chains, headgroups, and the solvent which results in a range of decay lengths spanning 0.8-2.5 8, and d,, values varying from 10 to 30 A. This situation is difficult and, in instances, impossible to correlate with some of these proposed explanations. However, the mechanism based on a thermally mobile and responsive bilayer is easily consistent with the measured behavior. The model system presented here clearly predicts a repulsive force of entropic origin over a wide range of parameter values. The separation dependence for the repulsive region in the two cases K , = Kh = 0.0546 J/m2, r = 0.0 and K, = 0.0546 J/m2, Kh= 0.437 J/m2, rq = 5.0 A is %own in Figure 7. Both sets of data have been fitted to an exponential with a decay length of roughly 2.1 A. Between b = 17.5 and 20 A, the decay is much faster than exponential, and at 20 A the pressure is definitely attractive. As no account of van der Waals potentials is considered here, this attraction stems from the two above-mentioned attractive sources: internal harmonic correlations between the particles as they cross the midplane49s50and charge correlations.'IJ2 Unfortunately, the form of the more rapid decay in the repulsion is difficult to follow by simulation as in this region the simulation error is relatively high. However, beyond 20 A we find, again with reasonable accuracy, a definite attractive force. Even in the absence of van der Waals forces, this region (17.5-20.0 A) marks an upper bound to the "swelling distance", d,,,,,. It is apparent from the figure that, to within simulation error, the above two cases give the same
Simulation of Flexible Zwitterions 2.0
E
1.5
P;,
1.0
2 5 Y)
&
0.5
8 4
0.0 -0.5
The Journal of Physical Chemistry, Vol. 95, No. 16, 1991 6359
1
i1
Q
cl
-1.0 -1.5 -0.5
I
I
I
5
10
15
20
Separation (A) Figure 7. Separation dependence of the pressure for the two cases Kc = Kh = 0.0546 J/m2, r = 0.0 (filled squares and solid line) and K, = 0.0546 J/m2, Kh = 0337 J/mZ, rq = 5.0 A (open diamonds and dotdashed line); the open triangles represent the uncharged case of K, = Kh = 0.0546 J/m2 and rq = 0.0. The symbols represent simulation points, and lines represent linear least-squares fit to the data. The particles have equal 2-A hard-core radii. The pressure is plotted on a density scale in molar units. The conversion factor to dyn/cm2 is 2.478 X IO7. TABLE III: Dependence on K, for fixed KI, = 0.437 J/m3 and rq = 5.0 A d KO = 0.0344 K, = 0.0546 K, = 0.0944
Po. M Po,dyn/m2 A,
A
d,,, A
J/m2 3.01 x 103
7.46 X 1Olo 2.3 22.5-25.0
J/m2 4.13 x 103 1.02 X IO" 2.1 17.5-20.0
J/m2 8.16 x 103
2.02 X loll 1.8 15.0-17.5
'A list of the least-squares fitting parameters, Po and A, of the exponential form Pi,,(6) = PO exp(-b/A). PO is quoted in both molar units and in dyn/cm2. The last row gives estimates of d,,, the equilibrium separation, quoted as a small range in separation: the lower (upper) end of the range is the largest (smallest) separation at which a repulsion (attraction) was reliably calculated. The particles have equal 2.0-A hard-core radii.
repulsive force and "swelling distance", suggesting some insensitivity to minor changes in the internal zwitterion structure, that is, with rq, Presumably, for widely different values of r9 divergences will appear, but such variations in rq are not consistent with a model of phospholipids. For the case K, = Kh = 0.0546 J/m2, r = 0.0, we include the force for uncharged particles 1 and 2. I n this region the calculated points are similar to the charged case. In fact, at separations less than 12.5 A the points coincide. At larger separations, however, the uncharged case shows a stronger repulsion. These observations lead us to conclude that between 12 and 20 A the force is dominated by restrictions in the movements of headgroups on one layer by the presence of the opposing layer. At still larger separations the attractive component caused by interlayer electrostatic correlations becomes significant and reduces the repulsion compared to the uncharged case. The log-linear fit to the data for the uncharged case gives a decay length of 2.6 A, i.e., slightly slower than for charged particles. The general conclusion of a survey of the effect of chain melting and chain heterogeneity on the extent of hydration is that increases in either temperature or degree of polyunsaturation in the chains result in increases in the hydration of bilayers, as measured by X and d,,,. Such changes in the state of the hydrocarbon core are reflected in values of the parameter K,: large values correspond to chains in the gel phase while smaller values imply a melted state. Keeping the average charge-charge distance within each zwitterion constant (with K,,= 0.437 J/m2 and rOs = 5.0 A, giving an average separation of 5.1 A), we can test the effect of varying the chain spring constant. Figure 8 shows the separation dependence of the pressure for three values of K,. In each of these cases an exponential force law can be 'fitted" to the data with decay lengths of 1.8, 2.1, and 2.3 A, respectively, for the K, values of 0.0944, 0.0546, and 0.0344 J/m2. The value of d,,, also increases proportionally (see Table 111).
5 5
10
15
20
25
Separation (A) Figure 8. Se ration dependenceof the pressure for the case of fixed Kh = 0.437 J/m p"and rps = 5.0 A for three values of K,: 0.0344 J/m2 (filled triangles and solid line), 0.0546 J/m2 (open diamonds and dot-dashed line), 0.0944 J/m2 (filled inverted triangles and dashed lines). The symbols represent simulation points, and lines represent linear leastsquares fit to the data. The prticles have equal 2-A hard-core radii. The pressure is plotted on a density scale in molar units. The conversion factor to dyn/cm2 is 2.478 X 10'. Table I11 summarim the least-squares fitting parameters for these cases. The open circles represent repulsive pressure values found by simulation of a case resembling the more rigid zwitterion system of ref 1 1 (see ref 63). TABLE I V Dependence on Hard-core Radii for Fixed K, = K, = 0.0546 J/m2 and rq = 0.0 & 2.0/2.0-A radii 2.0/3.0-A radii 1.4 x 104 Po, M 4.13 x 103 PO,dyn/m2 1.02 x 10" 3.46 X IO" A, A 2.1 2.5 dmx, A 17.5-20.0 25.0-27.5 OA list of the least-squares fitting parameters, Po and A, of the exponential form Pi&) = po exp(-b/A). po is quoted in both molar units and in dyn/cm2. The last row gives estimates of d-, the q u i librium separation, quoted as a small range in separation: the lower (upper) end of the range is the largest (smallest) separation at which a
repulsion (attraction) was reliably calculated. What must be true of both these simulations and of experiment is that the fitted exponential curves (the lines of Figure 8) provide only a rough average description of the forces over a particular separation range. This is a regime in which the loss of entropy due to the restriction of molecular fluctuations is dominant (true even in the absence of particle size), more or less as in the manner hypothesized by Israelachvili and Wennerstrom.'* One can rationalize this contribution further by considering the case of infinite KO Le., the rigid case with no perpendicular molecular fluctuations; a repulsion here is then due to the restriction of the orientational freedom of the zwitterions as the two layers overlap. Thus, necessarily, with a finite K, value there is both a "center-of-mass" freedom and an orientational freedom which are increasingly restricted as the two layers approach. Clearly, the perpendicular freedom allowed to the negative ions serves to extend the range of the repulsion! Indeed, this is not totally unlike the mechanism responsible for polymeric stabilization of colloidal suspensions,M except on a much smaller length scale. One must not forget that, for a given value of K,,at sufficiently large separations where the particle motions are negligibly restricted, attractive electrostatic correlations (and of course van der Waals forces) dominate. At smaller separations hard core effects are more apparent (although they do influence the force over all separations at which steric interaction occurs) as is clear from Figure 8, which shows the merging of the three curves, irrespective of the hydrocarbon "chain" (K,). Although the predictions of the model, as values of X and d,, summarized in Table 111, are consistent with measured trends under different chain conditions there is a current inadequacy in the theory. Different phospholipid chains have different (in(64) de Gennes,
P.Ado. Colloid Interface Sci. 1987, 27,
189.
6360 The Journal of Physical Chemistry, Vol. 95, No. 16, 1991 2.0 L
2
1.5
6
e
lI_
v1
0.5
E
3
+I
0.0
-0.5 -1.0
.,
1
'
1
'
1
'
1
\
; J
\. l
's
\9 '
,
I
,
I
.'
L3
I
(65) Parsegian, V. A.; Fuller, N.; Rand, R. P. Proe. Nod. Acad. Sci. U.S.A. 1979, 76. 2750.
Granfeldt and Miklavic
Conclusions A simple model for a zwitterionic molecule is presented as an extension of earlier work," now including a degree of motional freedom normal to the average hydrocarbon/water interface having as its physical origin a partial diffusive ability of the hydrocarbon chains out of the hydrocarbon. Using MC simulations, we have investigated the behavior of this model as concerns the conformational properties and the interlayer interaction. Upon choosing suitable representative systems, we have tested the model against some spectroscopic studies of phospholipid systems. From the comparison with and examination of experimental data, we can firstly conclude that our results agree qualitatively with experiment and secondly that the interpreted near-parallel orientation of the PN dipole (with an angle between 15 and 30°) is as much a consequence of the flutuation of the entire phospholipid in the direction perpendicular to the plane of the bilayer as of intermolecular electrostatic interactions. The interaction between two opposing layers is shown to be repulsive under conditions where the thermal motion of the particles is restricted by the presence of the opposing surface, "swelling" the two opposing layers to a finite distance, d,,, at which it is then balanced by attractive electrostatic correlation forces. The net repulsion is found to have an exponential behavior with magnitude and decay similar to fitted experimental values, and by varying suitable model parameters the simulations have shown that the model is qualitatively consistent with observed trends. The model of Israelachvili and Wennerstrbmls in which so-called hydration forces are explained as having an entropic origin is basically supported by the present calculations: the increase in free energy of interaction (loss of entropy) arises as a result of the restriction of both the center-of-mass freedom and orientational freedom of the molecules. As a general conclusion, the thermal agitation of amphiphilic molecules, in and out of hydrocarbon cores, provides a simple mechanism to account for the repulsions measured in such systems while remaining consistent with spectroscopic results of headgroup conformation. The nature of the headgroup modifies this underlying repulsion. For example, in our case we studied only particle size dependence, while in the case of the nonionic amphiphile, pentaoxyethylene dodecyl ether, the latently hydrophobic headgroup has an opposite effect, introducing an attractive contribution. Acknowledgment. We are extremely grateful to Dr. Bo J6nsson for his enthusiastic interest in this work and for his suggestions of improvement to the manuscript. We would also like to express our gratitude to Professor H. Wennerstrom and to the authors of ref 19 for providing a preprint of refs 18 and 19, respectively, and for subsequent discussions of their models and results.
