7201
J. Phys. Chem. 1994,98, 7201-7210
Monte Carlo Simulation of Lateral Distribution of Molecules in a Two-Component Lipid Membrane. Effect of Long-Range Repulsive Interactions Istvan P. Sugar,' Daxin Tang, and Parkson Lee-Gau Chong Departments of Biomathematical Sciences and Physiology t Biophysics, The Mount Sinai Medical Center, New York,New York 10029, and Department of Biochemistry, Temple University School of Medicine, Philadelphia, Pennsylvania 191 40 Received: January 31, 1994; In Final Form: May 11, 1994'
Recently a series of dips has been observed in the plot of E / M (the ratio of excimer to monomer fluorescence intensity) vs the mole fraction of l-palmitoyl-2-( l0-pyrenyl)decanoyl-sn-glycero-3-phosphoholine(Pyr-PC) in Pyr-PC/DMPC binary mixtures in the liquid crystalline phase at 30 OC (Tang and Chong Biophys. J. 1992, 63,903-910). In order to reveal connections between theE/Mdips and the lateral distribution of thecomponents, the membrane is simulated by means of Monte Carlo methods. The simulations confirm the assumption that pyrene-labeled acyl chains tend to distribute regularly in the membrane a t mole fractions where E / M dips are detected. The simulations also show that through the assumption of long-range repulsive interaction between the pyrene-labeled acyl chains (a) the area of regular arrangements have maxima a t the critical concentrations if the repulsion is strong enough, (b) in the case of weak repulsion the area of regular arrangements is a monotone increasing function of the concentration of the labeled chains, (c) when an E / M dip is measured a t critical concentration Xi, then a cluster of type i regular arrangement of labeled chains is always percolated throughout the membrane, (d) when an E / M kink is measured at Xi, a cluster of type i is rarely percolated, (e) between one critical concentration Xi and the next one &+I, an (order of type i ) (disorder) and then a (disorder) (order of type i+l) phase transition takes place.
-
-
Introduction A number of membrane activities, such as membrane permeability, membrane fusion, protein insertion, interbilayer lipid transfer, lipid flip-flop, and the activities of membrane proteins, have been linked to lipid lateral organization. The lateral distribution of the membrane components depends on external conditions such as temperature and pressure, on membrane composition, and on theinteractionsbetween the various molecular constituents. In the case of two-component bilayers there are three major types of lateral distribution of the constituents: domainsegregated, random, and regular. These types of distribution are strongly related to intermolecular interactions. In the case of domain segregated distribution the minor component tends to form clusters when similar components attract each other more strongly than do those that differ. These types of lipid membranes, thermodynamically classified as nonideal mixtures possessing positive deviation from ideality, have been studied inten~ively.3.~+'J7 Distribution is random when the interaction between similar components agrees with the interaction between components that differ. These systems are called ideal binary mixtures. There are several examples of almost ideal mixtures of two-component membranes, for example, when the head groups are identical while the lengths of the respective acyl chains differ by one or two methylene units.3.9J7 Lipids can be regularly distributed when different components attract each other stronger than the similar ones. According to the thermodynamical classification, these nonideal mixtures possess negative deviation from ideality. Somerharju et al.15 reported the first evidence that in the liquid crystalline phase 1-palmitoyl-2-(1O-pyrenyl)decanoyl-sn-glycero-3-phosphocholine (Pyr-PC) is regularly distributed in lipid membranes. They showed that plots of the ratio of the excimer to monomer
* Address all correspondenceto this author at The Mount Sinai Medical
Center. Phone: (212)241-5851. *Abstract published in Advunce ACS Absrrucrs, July 1, 1994.
0022-3654/94/2098-7201$04.50/0
fluorescence intensity (E/* vs the mole fraction of Pyr-PC in egg yolk lecithin and in L-a-dipalmitoylphosphatidylcholine (DPPC) vesicles have several linear regions separated by kinks. A theory has been established by Virtanen et aLZ1to determine the critical concentrations at which the kinks may be observed, assuming that at critical concentrations Pyr-PC will distribute regularly and form hexagonal superlattices throughout the membrane. Recently,in addition to kinks, we have clearly observed a series of dips, in the E/Mvs Xcurve for Pyr-PC/DMPC mixtures over a wide range of mole fractions above the phase transition temperatureofeachcomponent.19Thedips werelocatedat critical concentrations, in agreement with the hexagonal superlattice model, originally proposed by Vitanen et a1.**and recently further developed by Tang and Ch0ng.1~The geometric interpretation of the critical concentrations was proven to be very successful, correctlypredictingthe position of more than 12 dips. Compared to the result obtained by Somerharju et al.,15 this result provides more convincing evidence that lipids in Pyr-PC/DMPC can be regularly distributed at critical mole fractions. Lipid lateral distribution at noncritical mole fractions, however, is virtually unexplored; furthermore, from a thermodynamic viewpoint, a perfect hexagonal superlattice of Pyr-PC may form only at absolute zero temperature;this implies that, at room temperatures, regular regions must coexist with irregular regions. The relative contribution of these two distributions has not yet been studied. In this paper, lipid lateral distribution in Pyr-PC/DMPC twocomponent membranes is simulated at different mixing ratios, and at different intermolecular interactions by means of Monte Carlo (MC) methods, in order to reveal the real nature of the lateral distribution at raom temperature and at critical and noncritical Pyr-PC concentrations. Model
Virtanen et a1.21have developed a geometrical model for regular distributions of lipid molecules in two-component membranes. 0 1994 American Chemical Society
Sugar et al.
7202 The Journal of Physical Chemistry, Vol. 98, No. 29, 1994 a
O
Y
Y
a X
Y Y
Figure 1. Examples of perfect regular arrangements of pyrene-labeled chains. (a, Top left) Type 1 perfect regular arrangement at X = XI (= 0.666) critical Pyr-PC molar ratio. Lattice constants: u = 1, b = 2. (b, Top right) Type 2 arrangement at X = X2 (= 0.5). Lattice constants: u = 0, b = 2. (c, Bottom left) Type 3 arrangement atX = X3 (=0.286). Lattice constants: u = 1, b = 3. (d, Bottom right) Type 3 arrangement at X = X3 (=0.286). Lattice constants: u = 2, b = 3. Filled circles represent the labeled chains, while there are unlabeled chains at the other lattice points.
The model successfully predicts the critical concentrations of the Pyr-PC/ DPPC’5 and Pyr-PC/DMPC’g binary mixtures. We intend to develop a physical model of these two-component membrane systems by introducing interaction potentials between the molecules, keeping intact the geometrical characteristics of the Virtanen model. The Geometry of the Physical Model. The physical model simulates the equilibrium conformations of a monolayer of the bilayer membrane, assuming that the conformations of the monolayers are independent from each other. The assumption is based on an observation that the monolayers of a DPPC membrane are not thermodynamically c o ~ p l e d . ~Since ~ J ~ the acyl chains of membrane lipids have been shown to be situated on a triangular lattice,I2 it is assumed that the pyrene-labeled and unlabeled acyl chains of Pyr-PC/ DMPC two-component monolayers are also situated on a triangular lattice. The monolayer conformation, i.e. the distribution of the labeled acyl chains in the lattice, is represented by an N X N square matrix S (in this paper N = 12, 14, or 24). Each matrix element refers to a lattice point of the monolayer. In accordance with the lattice geometry, the following six matrix elements are the nearest neighbors of the 0th matrix element Sf,: S,-I,-I,SI-U,SI~I,S~+~~S~+VSI+~,+~. Each element of the S matrix has one of the two values SIj= 0 if the 0th lattice point is occupied by an unlabeled acyl chain, while Si, = 1 in the case of a labeled acyl chain. Thus the mole fraction of Pyr-PC is N
N
X=2cxSii/N2 is‘ j - I
Boundary Conditions. In order to reduce the edge effects in our small model system, periodic boundary conditionsZhave been implemented by assuming that replicas of the S matrix occur around the original S matrix. The abovedefinition of the boundary conditions makes it possible to calculate long-range interactions between lattice points even if the range is larger than the lattice
size.5 (Long-range repulsive interaction between the pyrenelabeled acyl chains is the physical basis of their regular distribution. This will be discussed at a later point in the paper.) Calculating, for example,the interaction between the ’ih lattice point and its mth right hand neighbor, we need to know the values of the respective S matrix elements. The state of the ijth lattice point is Si,, while the state of the mth neighbor is St.,,,odu+,,,a.The mod(k,l) function returns the remainder when the first argument is divided by the second argument. At zero remainder mod(k,l) = 1. Let us consider also a numerical example. In the case of a 12 X 12 lattice the state of the 22nd right hand neighbor of the (3,8)th lattice point is given by the following matrix element: S3,30 =: S3,mod(30,12)= S3,6; i.e. the 22nd neighbor is situated to the right in the 6th column of the second replica of the S matrix. Note that in most of the Monte Carlo simulationsof membranes only short-range van der Waals interactions between lipid molecules are modeled; i.e. nearest neighbor interactions are considered between the lattice points.’* Modeling long-range interactions is computationallymuch more demandingand always approximate.’ Critical Concentrations. In the case of an infinite triangular lattice, the labeled acyl chains can be arranged in a regular hexagonal order at certain mole fractions. According to Virtanen et al.,z’ these critical mole fractions of the Pyr-PC molecules can be calculated by means of the following equation:22 X(a,b) = 2/(a2 - ab
+ b2)
(2)
where a and b are nonnegative integers. For example, the four largest critical mole fractions of Pyr-PC are as follows: XI = X(1,2) 0.667,X2 = X(0,2) = 0.5, X3 = X(1,3) = X(2,3) = 0.286, X4 X(0,3) = 0.222. In Figure 1, perfect regular hexagonal arrangements of labeled chains are shown at different critical mole fractions. In this paper we refer to the regular arrangement of Xi critical mole fraction as “type i regular arrangement”. In the case of afinite triangular lattice with periodic boundary conditions, eq 2 has a limited validity. The pyrene-labeled acyl
The Journal of Physical Chemistry, Vol. 98, No. 29, 1994 7203
Lateral Distribution of Molecules
0.0 .00. 00.00 .00. 0.0
chainscan be arranged in hexagonal order only at a certain subset of the above critical mole fractions. For example, in the case of a 12 X 12 triangular lattice, regular hexagonal arrangements can be produced at XI,Xz, Xq, and X5 critical mole fractions but not at X3, while Xz and X , are the only critical mole fractions of a 14 X 14 triangular lattice. (One can try to build up a type 3 regular distribution at X3 critical concentration in the case of a 12 X 12 triangular lattice, to show that the order is always broken somewhere at the boundary between the lattice and its nearest replicas.) Thus if we want to simulate the monolayer conformations within the interval of two consecutive critical mole fractions, then the possible lattice sizes of the monolayer depend on the chosen interval. Conformational Energy. The energy of the monolayer at an S conformation is
5-2
s-1
0 0. 0000
. 0 . 0 .
O~OOO .
~O08O082 0 000
s-3
s=4
0.00 0000. .00000 000.000
where z = 6 is the coordination number. By means of the above equation one can eliminate every and Ms,from eq 3 and the conformational energy can be obtained as follows:
00.0 5-5
S=6
Figure 2. Different ranges of chain-chain interactions. The interaction of a chain with all the other chains is partitioned into interactions with chains in different ranges: e.g. s = 1, nearest neighbor interactions; s = 2, next nearest neighbor interactions; etc. In each figure the central chain (filledcircle)interacts with thesurroundingsixchains(fillcdcirclcs). Open circles represent other chains.
accepted if the following Metropolis criterion is fulfilled:
6E RAN < exp - -
kT
(7) is the cooperativity parameter of sth range interactions. The explicit expressions of E;/, and E:, interaction energies will be discussed after two subsections. Kawasaki Dynamics. Simulating the equilibrium fluctuations of the membranecauses theconformationalchanges to take place through elementary exchanges of nearest neighbor unlike chains. An elementary exchange is simulated by the followingsteps: (a) A pair of nearest neighbor unlike chains is selected randomly. (b) A trial membrane conformation is produced by exchanging the positions of the selected chains. (c) The trial conformation is
(8)
where RAN is a random number generated within a (0,l) interval, T i s absolute temperature, and k (= 1.987 cal/mol/deg) is the Boltzmann constant. From eq 6 it follows that the energy difference, bE, between the trial and original conformation is
Through the Metropolis criterion the membrane is driven toward equilibrium independently of the initial conformation, and eventually the energy of the fluctuating conformationsfollows a Boltzmann distribution. According to eq 9, the equilibrium conformations of the membrane depend on the cooperativity parameters W. Short- and Long-Range Interactions between Chains. In this section a cooperativity function w will be suggested by means of which one can calculate the cooperativity parameters w‘ at any s index. It is important to note that a characteristic distance of acyl chains R belong to each s index, and thus the cooperativity function w‘ can be given also in the function of these characteristic distances w(R). For example, if R is given in lattice units, then w ’ = w(l),wZ= w(&),w3=
where
0 00.
00000 000000 .00.00. 000000 00000 .00.
%oooo%oo~
In eq 3 only the conformation-dependent energy terms are considered. For example, the phosphatidylcholinephosphatidylcholine head grouphead group interactions are disregarded because their contribution to the monolayer energy is constant at any lattice conformation. The first two terms in eq 3 refer to the intrachain energies, and the other terms in the summation, to the interchain energies. NI and N , are the number of labeled and unlabeled chains, respectively,while Eland Enaretheenergies of a labeled and an unlabeled chain. In the calculation of the interchain energies the summation is taken over different ranges of interactions. For example, s = 1 refers to nearest neighbor interactions, while s = 2 to next nearest neighbor interactions. Figure 2 shows the interactions considered at different s values. In eq 3 XI,Ms,, and are the number of sth range interactions between labeled, unlabeled, and labeled-unlabeled chain pairs, and E!/, EL, and E:! are the respective interchain energies. In the case of periodic boundary conditions, the following relationships hold for every s index:*
00.0 .0000 00000. 000.000
w(2),w4=w5=w(fi),Hp=w(3)
(see Figure 2). The cooperativity function w is the sum of two terms, w , and~ wlw, referring to the short- and long-range chainchain interactions, respectively. The interaction between unlabeled hydrocarbon chains is inversely proportional to the 5th power of the interchain distance.1° This is a short-range attractive interaction (van der Waals interaction). The van der Waals interaction between labeled chains (or between unlike chains) follows the same power law, but the proportionality constant is different. Thus, the shortrange term of the cooperativity function is
-
where wshoflcan be calculated from the short-range interaction
7204
The Journal of Physical Chemistry, Vol. 98, No. 29, 1994
energies between nearest neighbor chains according to eq 7 and R is the interchain distance given in lattice units. The success of the geometrical interpretation of the critical concentrations of the Pyr-PC molecules1qshows us that a longrange repulsive potential between thepyrene-labeled hydrocarbon chains should be an essentialpart of the intermolecular interactions because (1) a repulsive interaction explains the tendency of the labeled chains to form hexagonal superlattices, and ( 2 ) a longrange interaction explains that the tendency of forming hexagonal superlattices persists at low Pyr-PC concentrations. Kinks (or dips) have been found in the experimental E / M vs X curves at X < 5%,15J9 where the average distance of the labeled chains is more than 6 lattice units (> 48 A).1z What could be the physical origin of this long-range repulsive interaction between the pyrene-labeled chains? According to the electronic structure in the ground state, pyrene is neutral and its permanent dipole is negligible. It is assumed, however, that the repulsive interaction is due to the presence of the bulky, perturbing pyrene moiety, which induces steric elastic strain in the alkyl chain 1atti~e.l~ Not knowing the form of the interaction potential, arising between the labeled acyl chains, in the MC simulations, we assumed a very simple form of the long-range term for the cooperativity function, as follows:
-
where wlongcan be calculated from the long-range interaction energies between nearest neighbor chains according to eq 7. In the case-of repulsive long-range interactions between - the labeled chains, wlongis negative. Throughout this paper (wlong( is referred to as "the strength of repulsion". Equation 1 1 defines a long-range interaction' because the interaction energy falls off no faster than R-d where d is the dimensionality of the system (in the case of membranes d = 2) and R is the interchain separation. Thus we choose the fastest decaying long-range interaction in our system. The technical advantage of this choice is that the infinite sum in eq 9, i.e. the energy difference,bE, between the trialand original conformationsofthe membrane, is finite. At long enough ranges, Ro(>>l), the summation in eq 9 can be replaced by the following integral:
where t is a unit vector between a pair of nearest neighbor unlike chains selected during an elementary step of the Kawasaki dynnamics; R is the vector between the selected labeled chain and a point of the membrane of polar coordinates R and cp. The origin of the polar coordinate system is defined by the selected labeled chain. p(R) is the local surface density of the labeled chains in a point of the membrane at R. In the Appendix it is pointed out that the absolute value of the integral in eq 12 is bounded by 47&lon#Ro. In the actual computations long-range interactions have been considered up to the 30th neighbor, i.e. Ro = 31. With this cutoff length the error of 6E is less than 3%. Note that the interaction energy between a selected labeled chain and the other labeled chains diverges logarithmically with the upper limit of the integration. However, this divergence does not affect the goodness of the Monte Carlo simulation because theenergydifference bE rather than theinteraction energymatters at every decision making. Monte Carlo Cycle. An MC cycle simulates the thermal fluctuations of the monolayer within a very short time interval. During the cycle, pairs of nearest neighbor chains are randomly selected z ( N / + N n ) / 2times, each time they get an opportunity to exchange position (where z(N1 + N n ) / 2is the total number of nearest neighbor pairs in a lattice). This way, on average, each pair of chains is selected once during an MC cycle.
Sugar et al. SMpPhot Analysis. At the end of each MC cycleseveral aspects of the membrane conformation have been analyzed. From the number of nearest neighbor labeled chains Ni, we could determine the following quantity, which is proportional to the measurable E / M ratio:
E/M
-
(N,)/N,
(13)
where the average of Nj, and Nl is proportional to the probability of excimer and excited monomer formation, respectively. In order to obtain information about the type and amount of regular distribution of the labeled chains, each snapshot has been transformed by means of a pattern recognition algorithm. When the molar ratio of the Pyr-PC moleculesXis between consecutive critical concentrations X,, and X,,+l, the pattern recognition algorithm searches for type n and type n+l regular hexagonal arrangements. For example, Figure 3a,b shows snapshots (i.e. structure matrices) taken at Pyr-PC molar ratios ofX = 0.52 and X = 0.6, respectively. These are mole fractions between critical concentrations XI(= 0.667) and X2(= 0.5). Diagonal crosses and points represent labeled and unlabeled chains, respectively. Investigating these snapshots, one can recognize islands where labeled chains form type 1 and type 2 regular hexagonal arrangements. This visual recognition of certain hexagonal arrangements is automatized by means of our pattern recognition algorithmwhich transforms thesnapshot (Le. theactual S matrix) to a pattern matrix P. The elements of the pattern matrix are defined as follows: Pi, = 1 if the ijth lattice point is within an island of type I hexagonally arranged labeled chains. Pij = 2 if the ijth lattice point is within an island of type 2 hexagonally arranged labeled chains. Pij = 0 if the ijth lattice point is out of the areas where the above two types of regular arrangements can be found. In our example Figure 3c,d shows the pattern matrices obtained from the snapshots in Figure 3a,b, respectively. In these pattern matrices stars and crosses represent areas where regular hexagonal arrangements of type I and type 2 can be found, while none of these two regular arrangements was found at the blank areas. The pattern matrix directly characterizes the types and amount of the regular distributions in a snapshot. The pattern matrix, P,is further analyzed by determining the number of lattice points belonging to an area of certain regularity, and the number of isolated areas (Le. clusters) of a certain regularity. This analysis of the pattern matrix is performed by means of a common cluster counting algorithm16 which enables us also to recognize large clusters which span the pattern matrix either horizontally or vertically. A pattern matrix with such a large cluster is called percolated (see e.g. the largest type1 cluster in Figure 3d).
Results and Discussion MC simulations of the Pyr-PC/DMPC membranes have been performed in the [XI&] and (Xz,X3] intervals of critical concentrations with lattice sizes of 12 X 12 and 14 X 14, respectively. (Square bracket and ordinary parenthesis mark the closed and open ends of an interval, respectively.) Simulations were of repulsion (Le. at - carried out at five different strengths lwlongl= 0-4kcal/mol/chain), while wshort= 0 was taken always. In every simulation, the temperature was set to 300 K. Longrange interchain interactions have been calculated up to the 30th neighbors of each labeled chain. Each simulation has been started from random distribution. The attainment of equilibrium was monitored by the lattice energy, and it was usually attained after 500 MC cycles. After equilibrium was attained, snapshots of 1000consecutiveMCcycles wereanalyzedby meansof the pattern recognition and cluster counting algorithms to obtain statistical averages characterizing the monolayer conformations. Checking the effect of lattice size, we also ran simulations on a lattice of 24 X 24, from X = 0.50 to 0.667 mole fraction at
Lateral Distribution of Molecules
The Journal of Physical Chemistry, Vol. 98, No. 29, 1994 7205 " X X .
A.
, .
X
+++*
**
C.
*it+ tt +t.*
tt**
..
,; +.* ++( * * .lt*. +
+ '
*
+
+ + e *
,++t + + * . **
..*..
*
**.
+*
++
.
+*
++
t.
i. t
++ *
...................
*.
*."*,+ +'t
. +; +.++. ;*;.;**
. I . .
++
** ++t . ** ** * *
' *******
........ I.
2 kcal/mol/chain, while the increasing number of type 1 clusters at = 2 kcal/mol/chain shows a slight dominance of nucleation over the coalescence. From X = 0.54 the coalescence of type 1 clusters becomes the leading phenomenon, and as the critical concentration XI is approached, a percolated cluster of type 1 appears with a high frequency (Figure 7b). Figure 6c shows that in the case of strongrepulsion the number of clusters of irregularly arranged pyrene-labeled chains has a maximum in the [X~,XI]interval. At different strength of repulsions the position of this maximum corresponds with the
&,,.el
Lateral Distribution of Molecules half-percolation point of type 1 clusters and - clusters of irregular arrangements as well. For example, at wlm8= -3 kcal/mol/ chain the maximum is at X = 0.55-0.58, while the percolation frequency of type 1 clusters is 0.5 at X = 0.56 (Figure 7b); for clustersof irregular arrangements the half-point of the percolation is at X = 0.58 (Figure 7c). Thus, increasing the concentration from X2 causes the coalescenceof small clusters of type 1 to disconnect the percolated cluster of irregular arrangements, resulting in an increase in the number of clusters of irregular arrangements. Above the percolation point of type 1 clusters, however, the shrinking and disruption of the irregularly arranged areas takes place; i.e. the number of clusters of irregular arrangements decreases. E/MCurves and Cluster Percolations. By comparing the above results of the pattern distribution analysis and thecurves in Figure 4, we find that in every case qualitative - changes take place when the repulsion becomes strong, i.e. lwlonJbecomes larger than 2 kcal/mol/chain. Now, is there any property of the pattern distribution connected to the E / M dip formation? We know from Figure 4 that at the same wlongone can get dip at X I critical concentration, but there is no dip at X2; why is this? These questions can be answered by examining the percolation frequencies of the clusters of different patterns (see Figure 7). Figure 7b shows that in the case of strong repulsion there is almost always a percolated cluster of type 1 arrangement near X I critical concentration (Le. the percolation frequency is very close to 1 ) . This extended cluster of regularly arranged labeled chains results in small excimer formation probability, Le. E / M dip. On the other hand, the clusters of type 2 arrangements are very rarely percolated at X2 critical concentration (see Figure 7a), and thus the excimer formation probability is not low enough to get an E / M dip. The percolation frequency of the clusters of irregularly arranged pyrene-labeled chains drops to almost zero simultaneously with the percolation of the type 1 clusters (see Figure 7c). Series of Concentration-InducedPhase Transitions. According to our Monte Carlo simulations, the formation of E / M dips is associated with a series of changes in the lateral distribution of the labeled chains in the membrane. Increasing the concentration from critical concentration Xi to Xi+l causes two qualitative changes to take place in the lateral distribution. Close to Xi, a cluster of type i regular arrangements of labeled chains is percolated throughout the membrane; however, as the concentration of the labeled chains increases, the percolation frequency suddenly drops to zero. On the other hand, as the Xi+l critical concentration is approached, the percolation frequency of type i+ 1 clusters suddenly increases from 0 to 1 . These sudden, but still continuous, changes in the percolation frequencies signify concentration-induced phase transitions in the membrane. The type of these phase transitionscannot be determinedaccording to the Ehrenfest’s classification, Le., by analyzingthe singularities of the functions of statistical averages, because these functions never have singularities in the case of small systemsa6 In the case of small systems the type of the phase transition can be determined from the distribution functions of the fluctuating extensive parameters of the system: in the case of first-order phase transition the distribution function is inhomogeneous at the transition region, while in the case of second-orderphase transition it is homogeneou~.~ In our Monte Carlo simulations at wlong= -3 kcal/mol/chain, the percolation frequency of type 1 regular arrangements changes suddenly in the concentration region of (0.55,0.6).In this region the distribution function of the relative area of type 1 regular arrangements was found to be homogeneous (not shown); i.e. a second-order phase transition takes place. The series of concentration-indud phase transitions, predicted from our Monte Carlo simulations, are supported by recent experimental EIMdata (Somerharju et al., University of Helsinki,
The Journal of Physical Chemistry, Vol. 98, No. 29, 1994 7209 Finland) showing enhanced fluctuations close to and on both sides of every critical concentration. The enhancement of the fluctuations at the phase transitions may explain our observation20 that dip positions are reproducible while the shape of the E / M vs X curve is less certain. Finally, we have to mention a two-dimensionalsystem, namely a submonolayer film, adsorbed on a periodic substrate, which also undergoes a seriesof concentration induced phase transitions. Within a certain temperature range, increasing the submonolayer density causes “commensurate” solid phases to alternate with “floating” solid phases.ll In a commensurate solid phase, an adsorbate lattice or superlattice is locked to the substrate periodicity. The phenomenological theory of these supported one-componentfilms’ may help us todevelop a phenomenological theory of our two-component membrane system, where the unlabeled chains (matrix lipid) correspondto the periodicsubstrate and the labeled chains to the submonolayer. Conclusions
In conclusion, the Monte Carlo simulation of the lateral distribution of molecules in a two-component membrane shows that, with the assumption of a long-range repulsive interaction between the labeled chains (eq l l ) , (a) between one critical concentration Xi and the next one X,+l an (order of type i ) (disorder) and then a (disorder) (order of type i + l ) phase transition takes place, (b) the areas of regular arrangements have maxima at the critical concentrations if the repulsion is strong enough, (c) in the case of weak repulsion the area of regular arrangements is a monotone increasing function of the concentration of the labeled chains, (d) when an E / M dip is measured at critical concentration Xi, then a cluster of type i arrangement of labeled chains is always percolated throughout the membrane, (e) when an E / M kink is measured at Xi, a cluster of type i is rarely percolated.
-
-
Acknowledgment. The authors thank Dr. J. Virtanen for valuable discussions and Ms. T. Hill for her contributions to the text. This research was supported by the U S . Army Research Office. A generous grant of computer time is acknowledged from the University Computer Center of the City University of New York. This research was conducted using the resources of the Cornel1 Theory Center, which receives major funding from the National Science Foundation and the IBM Corp., with additional support from New York State Science and Technology Foundation and members of the Corporate Research Institute. Appendix
In this appendix we estimate the integral in eq 12
where t is a unit vector between a pair of nearest neighbor unlike chains selected during an elementary step of the Kawasaki dynamics; X / 2 is the mole fraction of the labeled chains; R is the
7210 The Journal of Physical Chemistry, Vol. 98, No. 29, 1994
vector between the selected labeled chain and a point of the membrane of polar coordinates Rand q. The origin of the polar coordinate system is defined by the selected labeled chain. p(R) is the local surface density of the labeled chains in a point of the membrane at R. According to the above derivation, the up_per bound of the absolute value of the integral in eq 12 is 4rlw,,,l/Ro. References and Notes (1) Allen, M. P.;Tildesley, D. J. CompurerSimulutionofLiquids; Oxford Science Publications: Oxford, 1990. (2) Born, M.; Von Karman, Th. h e r Schwingungen in Raumgittern. Phys. Z . 1912, 13,297-309. (3) Brumbaugh. E. E.; Huang, C. Parameter Estimation in Binary Mixtures of Phospholipids. Methods Enzymol. 1992, 210, 521-539. (4) Von Dreele, P. H. Kstimation of Lateral Species Separation from Phase Transitionsin Nonideal Two-Dimensional Lipid Mixtures. Biochemistry
1918,17, 3939-3943. (5) Ewald, P. Die Berechnung Optischer und Elektrostatischer Gitterpotentiale. Ann. Phys. 1921, 64, 253-287. (6) HiI1,T. L. ThermodinumicsofSmulISyste~; Benjamin: New York, 1963. (7) Hill, T. L. Cooperutivity Theory in Biochemistry; Springer-Verlag: New York, 1985; Section 10. (8) Huang, K. Sturisticul Mechanics; John Wiley and Sons: New York, 1963. (9) Mabrey, S.; Sturtevant, J. M. Investigation of Phase Transitions of
Lipids and Lipid Mixtures by High Sensitiity DSC. Proc. Nurl. Acud. Sci. U.S.A. 1916, 73, 3862-3866. (10) Nagle, J. F.; Wilkinson, D.A. Lecithin Bilayers. Density Measurements and Molecular Interactions. Biophys. J. 1978, 23, 159-175. (1 1) Nelson, D. R.; Halperin, B. I. Dislocation-Mediated Melting in Two Dimensions. Phys. Rev. B 1979, 19, 2457-2484.
Sugar et al. (12) Ruocco, M. J.; Shipley. G. G. Characterization of the Subtransition of Hydrated Dipalmitoyl Phosphatidylcholine Bilayers. Kinetics, Hydration and Structural Studies. Biochlm. Blophys. Acta 1982,691, 309-320. ( 13) Schmidt, C. F.; Barenholz, Y .;Huang, C.; Thompson, T. E. Monolayer Coupling in Sphingomyelin Bilayer System. Nurure 1978, 271, 775-777. (14) Sillcrud, L. 0.; Bamett, R. E. Lack of Transbilayer Coupling in Phase Transitions of Phosphatidylcholine Vesicles. Biochemistry 1982, 21, 1759-1776. (15) Somaharju,P.J.;Virtamn,J.A.;Eklund,K.K.;Vainio,P.;Kinnunen.
P. K. J. l-Palmitoy1-2-pyrenddccanoylGlyccrophospholipids as Membrane Probes: Evidence for Regular Distribution in Liquid Crystalline Phosphatidylcholine Bilayers. Biochemistry 1985, 24, 2773-278 1. (16) Stauffer, D. Introduction to Percolution Theory; Taylor and Francis: London-Philadelphia, 1985. (17) Sugar, I. P.;Monticelli, G. Interrelationships between the Phase Diagrams of the Two-Component Phospholipid Bilayers. Biophys. J . 1985, 48,283-288. (18) Sugar, I. P.; Biltonen, R. L.; Mitchard, N. Monte Carlo Simulations
of Membranes: Phase Transition of Small Unilamellar Dipalmitoylphosphatidylcholine Vesicles. Methods Enzymol. 1994, 240, 569-593. (19) Tang, D.; Chong, P. L.-G. E/M Dips. Evidence for Lipids Regularly Distributed into Hexagonal Super-Lattices in Pyrene-PC/DMPC Binary Mixtures at SpcCific Concentrations. Biophys. J . 1992, 63,903-910. (20) Chong, P. L.-G.; Tang, D.; Sugar, 1. P. Exploration of Physical Principles Underlying Lipid Regular Distribution: Effects of Pressure, Temperature and Radius of Curvature on E/M Dips in Pyrene-Labeled FC/ DMPC Binary Mixtures. Biophys. J. 1994,66, 2029-2038. (21) Virtanen, J. A.; Somerharju, P.; Kinnunen, P. K. J. Prediction of Patterns for the Regular Distribution of Soluted Guest Molecules in Liquid Crystalline Phospholipid Membranes. J. Mol. Electron. 1988, 4, 233-236. (22) In contrast to the equation derived by Virtanen et a1.,2' eq 2 has a negativetermin thcdenominatorbccausehere thecoordinatesystemattached to the lattice is defined differently; i.e. for convenience the direction of the Y-axis has been reversed.
