- Email: [email protected]
Lateral diffusion in a binary lipid system by a computer simulation model
PHYSICA ELSEVIER
Physica A 223 (1996) 309-320
Lateral diffusion in a binary lipid system by a computer simulation model R.B.
Pandey
Department of Physics and Astronomy, The Program in Scientific Computing, University of Southern Mississippi, Hattiesburg, MS 39406-5046, USA
Received 4 October 1995; revised 19 October 1995
Abstract A computer simulation model is used to study the lateral diffusion in a two-dimensional binary mixture of large and small particles. The small particle represents the phospholipid while the large particle (square) models the planar macromolecule in the binary mixture. The variation of the rms displacement with time of these molecular species are studied as a function of their concentration and their diffusion constants are examined. The concentration of the large particles (pl), small particles (p2), and the vacancy (P3) are varied at the relative area fraction Ae = p3/(P2 q- P3) = 0.30 (low) and 0.60 (high) free area fractions. The diffusion constant of the small particles decays on decreasing their concentration - this decay signals the possibility of a transition between a fast diffusive regime (at high concentration) and a slow diffusive regime (at low concentration p2). The diffusion constant of the large molecules on the other hand shows an opposite behavior, it decays linearly on increasing their concentration (pl). The diffusion constant of the large molecule decays with their size while that of the small particles remains less affected by varying the size of the large molecules. The decay of the diffusion constant of the large molecule with its linear size is found to be nonlinear.
1. I n t r o d u c t i o n The molecular mobility of phospholipids in cell membranes plays a very important role in governing their physical and chemical properties [ 1 ]. Of our particular interest is the lateral diffusion of lipids in biphasic phospholipid monolayers [2,3]. The system consists of a binary mixture of two species - small phospholipids and large lipid macromolecules with their weight ratios of the order of 50:1 to 100:1. Extensive experimental studies [3,4] such as surface pressure measurements with FRAP method and fluorescence microscopy have been recently carried out to measure the diffusion constant of such binary mixtures as DLPC (L-alpha-Dilaroylphosphatidylcholine). Several 0378-4371/96/$15.00 (~) 1996 Elsevier Science B.V. All rights reserved SSDI 0378-4371 (95)00366-5
310
R.B. Pandey/Physica A 223 (1996) 309-320
interesting findings have been reported on the variation of the diffusion constant with the concentration of the constituents, their size, and the nature of the interfacial substrate among other parameters. The mobility of the two species obviously depends on several factors such as concentration of these molecules, free surface/volume, their relative size, shape, interaction, and temperature, etc. Understanding the diffusion in such a complex mixture is one of the most difficult theoretical problems due to strong correlations from small to large scales. Even the approximate theoretical methods such as mean field theories including selfconsistent field and effective medium [ 5] approaches are limited to monodisperse systems; most of these methods are known to be inadequate for taking into account the effects of correlations. Therefore, computer simulation methods remain one of the very important tool to understand the diffusion in such complex systems. Computer simulation models have been recently used to study the phase separation and phase transition in simplified amphiphilic systems [6,7]. The diffusion of amphiphiles are particularly used to monitor the evolution of micelle formation as the system evolves through a bicontinuous phase in a model of micoemulsion [7]. Because of the limitations on the computing resources (memory and speed), it is not feasible to take into account all the relevant details of such a binary system. Therefore, we will limit ourselves to simplified model. We study here the diffusion in a binary lipid system consisting of planar macromolecules and a small phospholipid by a simplified computer simulation model.
2. Model
We consider a discrete L x L square lattice. The small phospholipid molecule is modelled by a point-like particle while the large lipid macromolecule by a solid square of size l x I. Thus, a small particle occupies a lattice site and a large particle 12 lattice sites. We would like to initialize the lattice by randomly distributing these particles with a concentration Pl of large, and P2 of small particles; the concentration P3 of vacancy (i.e. empty sites) is P3 = 1 - p l - P 2 . A technical problem regarding the random packing arises when we attempt to distribute the large particles beyond a certain concentration (Plc), known as the first jamming/saturated coverage [8,9] or percolation threshold [ 10]; generally Plc depends on the size and shape of the objects, and decreases with the size of the object with power-laws for several well-defined shapes [ 8,11 ]. However, in order to get around this problem, one may distribute the large particles in a regular fashion, and then move them around randomly for a predetermined period of time before quenching them to obtain their random distribution. Such process of annealing followed by quenching may not provide a complete random distribution, but offers an alternate solution to jamming coverage. Our system is close to an annealed system (see below), therefore, we believe that this method of initializing our lattice system does not affect our results. Thus, we first distribute large particles on a fraction Pt of the lattice sites and then a fraction P2 of the lattice sites are randomly occupied by small particles. An excluded volume (hardcore) interaction is considered between the particles; a large
R.B. Pandey /Physica A 223 (1996)309-320
311
particle occupies l 2 lattice sites while a small particle occupies a lattice site. A lattice site cannot be occupied by more than one particle at a time. We stir the system by moving these particles randomly for a period of time as follows. For small particles, we choose a particle and one of its neighboring site randomly; if the neighboring site is empty then we move the particle. However, if the neighboring site is occupied by a small particle or by part of a large particle, then the particle remains at its old position. Essentially the same procedure is applied in moving a large particle as well. We attempt to move all the sites of an entire macromolecule by one lattice constant in a randomly chosen direction. Thus, for a successful move, there must be at least l neighboring sites vacant in the direction of the move. The probability of finding l empty neighboring sites along the edge of the large particle in the direction of the hopping is usually small in a dense system. The motion of the large particles is therefore expected to be slow. Furthermore, even though the center of mass of the large particle is shifted by one lattice constant in a successful hopping, l sites simultaneously become vacant and occupied in moving a large particle. Therefore, we attempt to move the large particle with a hopping rate ~ 1 / l (see below). The attempt to move each particle regardless of its success to hop to its neighboring site/sites is considered in defining the time step. The definition of Monte Carlo Step is generally arbitrary [ 12] ; an attempt to move each particle once is usually defined as one Monte Carlo Step (mcs). In our case here, there are two particles of different sizes, therefore we define a unit time step in which each small particle attempts to move once while a large particle moves with a probability 1 / l . A periodic boundary condition is used to move the particles across the boundaries. We consider two systems. ( 1 ) A quenched system, in which we move all the particles in the system for a pre-assigned period of annealing time, and then we quench the large particles. The "so-called" random distribution of large particles thus achieved in such an annealing process acts as a quenched random barrier for the mobility of the small particles. Such quenched porous media may be compared to a percolating porous media. However, one should be careful in comparing the physical properties of such porous medium with that of ordinary site/bond percolation, since not only the percolation threshold of such systems are different but also the morphology, and, therefore, the geometrical correlations may differ substantially. (2) An annealed system, in which we let the particles move around for a sufficiently long time so that we achieve their asymptotic transport behavior. The whole simulation is repeated for a number (N~un) of independent samples to obtain a reliable average of the transport quantities. During the simulation we keep track of the displacement of all particles and their center of mass, and we evaluate the rms displacement of each particle periodically during the entire course of time. The leading power-law behavior of the rms displacement (R) of a particle with time t is described by R = 2Dt k ,
(1)
where k is a power-law exponent, and 2D is a prefactor. For most annealed systems, k = 1/2 in the asymptotic regime, and the stochastic motion of the particles is called
R.B. PandeylPhysica A 223 (1996)309-320
312
,,,,/
(a)
10098" 76-
E
s-
~ 131
4-
~
3-
8
10-
| I
I
i
1 04
1 0s Time
I
(b) 2.22 -
2.20 -
aI
c
°
2.18-
A All ~
ffJ
& j
2.15-
All
n-
•
Vo°
•.
V..,~ .V~ V ...V
AI !re ~1 "~o
2.14-
I
* •
•
Vo
V o
~o 2.10-
6x 10 4
1 05 Time
F i g . 1.
R.B. Pandey/Physica A 223 (1996) 309-320
313
Fickian or diffusive with the Einstein diffusion constant D. The motion o f a particle can, however, be anomalous and nondiffusive in many systems such as in random percolating system at percolation threshold [10,13,14] or in a complex mixture o f interacting systems [ 15] where the exponent k v~ 1/2. In the latter case, the "diffusion constant" terminology for the prefactor ( D ) may be misleading. We will, however, stick to the diffusion constant [16,17] since the motion is diffusive ( k = 1 / 2 ) for most o f the concentration regime considered here. The diffusion constant D depends on the heterogeneity o f the host space [ 14] and the concentration o f the mobile species [ 16] which is the main focus o f our study in the following. We w o u l d like to mention that there are different methods o f estimating the diffusion constant. In continuum simulations, particularly in Molecular Dynamics, one usually evaluates the velocity autocorrelation function, and averaging over time then gives an estimate o f diffusion constant [ 18]. Alternatively, since the motion is diffusive, the instantaneous derivative o f the average value o f R 2 with time will give us the diffusion constant which should approach a constant value in the asymptotic limit ( t ~ cx~). The intercept o f a l o g - l o g plot o f the rms displacement versus time may also give us an estimate o f the diffusion constant which may vary somewhat with the number o f data points particularly when the data are fluctuating. In order to avoid confusion, we will concentrate to an effective diffusion constant ( 2 D ) , which is obtained by evaluating R 2 / t , where t is the largest time in our simulations. Since we are interested mainly in the variation o f the diffusion constant with the concentration o f the constituents rather than its absolute value, the evaluation methods should not affect our findings. Therefore, we will use the simple method here. Furthermore in order to reduce the relative fluctuations and to maintain the consistency, the period o f simulation and time t at which the diffusion constant is evaluated is kept the same for all the simulations to study its relative variation with the concentrations o f the two types o f particles.
3. Results and discussion The simulations are performed on square lattices o f size 128 x 128 and 256 × 256 with various sizes o f the large particles, 82 to 322 . The large particles are used to model the effects o f large macromolecules and their m a s s / s i z e reflects a relative measure with respect to small particles which represent the small phospholipids. The number o f independent runs Nrun are limited and varied with Nrun = 5 - 2 0 in order to optimize the computing resources and to reduce the statistical fluctuations. As we mentioned
Fig. 1. RMS displacement of small particles versus time on a log-log scale with (a) on the full time scale while (b) is only in the long time regime. Lattice of size 128 x 128 was used with the large particles of size 82 with 20 independent samples. The relative free area fraction Ae 0.30 at the concentrations Pl = 0.1,p2 = 0.63 (circle), Pl = 0.2,p2 = 0.56 (square), Pl = 0.3,p2 = 0.49 (triangle), Pl = 0.4,p2 = 0.42 (inverted triangle), Pl = 0.5,p2 = 0.35 (parallelogram). Open symbols represent the rms displacement of the particles in quenched system while the filled symbols represent the rms displacement in the annealed system. =
R.B. Pandey/Physica A 223 (1996) 309-320
314
0
1.00
0.98 o 1= n
0
0.96 -
0,94 -
8 0.92 a
) )
0.90
o Z
O 0~88 •
0.0
I 0.1
I 0.2 Concentration
I 0.3
• 0
of Small
0.4
I 0.5
I 0.6
Particles
Fig. 2. Normalized diffusion constant of small particle versus their concentration (P2) at the free area fraction Ae = 0.30. The statistics is the same as in Fig. 1. The open circle represent the diffusion constant of the particle in the quenched system while the filled circle, in the annealed system. The diffusion constants are normalized by their maximum value.
before, we consider two systems for comparison. (1) The quenched system in which the large macromolecules are quenched and act like quenched barriers for the mobile phospholipids. (2) The annealed system, in which both large and small molecules are mobile with hopping attempt rates inversely proportional to their linear size. We limit ourselves to two values of the free area fraction, Ae = P2/(Pz + P3), a low free area (Ae = 0.30) and a relatively high free area (Ae = 0.60). For both sets of free area fractions, the concentration of the large macromolecule (Pl), small phospholipids (P2), and the free area (P3) are varied with Pl + P 2 q-P3 = 1. All the data are generated on a Cray YMP and a Cray J916 machines; about 100 hours of CPU are used. We have studied the variation of the displacement and the root mean square (rms) displacements of both small and large particles and their center of mass. We monitor the displacements primarily as an alternate method to check whether the system has reached the asymptotic regime. The variation of the rms displacement with time is the main quantity we analyse here. The large particles move relatively slow, therefore we will mainly concentrate on the rms displacement of the small particles except toward the end to discuss the effects of the size of the large particles on their transport. The variation of the rms displacement of the small particle with time, at a free area fraction, Ae = 0 . 3 0 , is presented in Fig. 1 for both quenched and annealed systems for various concentrations of the constituents. From a first glance on Fig. la, it appears that there is no effect on
R.B. Pandey/Physica A 223 (1996) 309-320 1.00
315
-
•
0
©
g
© 0.95
-
II.
0
t~
E
0.90-
"6 c
0.850
i5 0 . 8 0 lo
g
o 0.75 z
0.70
0.00
I
I
I
0.05
0.10
0.15
Concentration
• I
I
I
I
0.20
0.25
0.30
0.35
of S m a l l
Particles
Fig. 3. Normalized diffusion constant versus concentration of the small particles, similar to Fig. 2, for the free area fraction Ae = 0.60. The statistics is the same as in Fig. 2.
changing the concentration ( p 2 ) of the small particles with a little difference in the temporal variations of the rms displacement for the quenched and the annealed systems. All the data seem to collapse nearly to a single curve with a linear diffusive power-law fit, i.e. R N t l / 2 especially in the long time (asymptotic) regime; from Eq. (1), k = 1/2. Even the intercept of these data points look nearly the same in Fig. la. However, a close examination of these curves (Fig. lb) reveal the fact that the diffusion constant (i.e. the intercept) depends on the concentration of the constituents; diffusion constants at some concentrations in the annealed system differ from those in the quenched systems. As we mentioned before, there are several methods to evaluate the diffusion constant. We estimate an effective diffusion constant De = 2D = ( R 2 ) / t as t ---+~ , the largest time in our study. In Fig. 2 we present the variation of the effective diffusion constant of the small particles with their concentration at the low free area fraction (Ae = 0.30). As expected, the diffusion constant for the small particles in the annealed system has comparatively less fluctuations than the corresponding values in the quenched system. In the annealed system, the diffusion constant decays on reducing the concentration of the small particles until a characteristic concentration (P2c) below which it remain nearly a constant to a low value. This may be a signal of a transition, from a relatively fast diffusion at/)2 > P2c to a slow diffusion at p2 < p2c. This transition may be different from the conductor-insulator transition in a static percolation [7]. The former is a dynamic heterogeneous system while the latter is a quenched inhomogeneous
R.B. Pandey/Physica A 223 (1996) 309-320
316 t.0-
0h91:: e5 ¢1 ¢D
O
0.8o c e~
0.7-
8
g 0.6
ca N
0.5
O
z
0.4
I
0.2
I
I
I
0.4 0.6 0.8 Concentration of Large Particles
I
1,0
Fig. 4. Normalized diffusion constant versus concentration of the large particles at the free area fraction Ae = 0.30 (open circle) and Ae = 0.60 (filled circle). The statistics is the same as the previous figures.
system. Unfortunately, a considerable amount of more computing resources are needed to quantify the possibility o f such a transition in our annealed system. The diffusion constant in our quenched system, on the other hand, shows a peak before falling off at nearly the Same value o f P2 as the diffusion constant i n the annealed system. However, since there is m o r e fluctuations in the quenched case, we should be Cautious in such conclusions for the quenched system. We have carried out a similar investigation also at a relatively large free area fraction (Ae = 0 . 6 0 ) . The variation o f the corresponding diffusion constant of the small particles with their concentration is shown in Fig. 3. It is relatively harder to evaluate the diffusion constant more accurately due to large fluctuations in the data here than that with low free area fraction. However, as with the low free area fraction, the data with the annealed system shows less fluctuations than that in the quenched system. We see a decay in the diffusion constant on reducing the concentration (P2); the peak at the concentration p2 "~ 0.30 may be an artifact o f fluctuations. The trend in the variation of the diffusion constant with the concentration in the quenched system is even less predictable which, nevertheless, seems to suggest an opposite trend to that in the annealed system. The diffusion o f large molecules is studied in the annealed system. The variation o f the diffusion constant with their concentration is presented in Fig. 4 at the free area fractions Ae = 0.30 and 0.60. Except at the low concentrations (pl = 0.1 - 0.4) at Ae = 0.30, the diffusion constant seems to decay linearly on increasing the concentration
R.B. Pandey/Physica A 223 (1996) 309-320
317
n
100E
9-
"6
B-
"E
~-
OU
7-
E
8
6-
== co
O
n:
[]
•
Q
0 1 04
1 05 Time
Fig. 5. RMS displacementof small particle versus time on a log-log scale at the free area fraction Ae = 0.30 (circle) and Ae = 0 . 6 0 ( s q u a r e ) . The size of the large particle is 82 in all the runs. Lattice of size 128 x 128 (open symbols) is used with 20 independentsamples while the lattice of size 256 x 256 (filled symbols) is used with 5 and 15 independentruns for Ae = 0.30 and 0.60, respectively. of large particles. Note that the diffusion constant (D1) of the large particle shows an opposite trend to that for the small particle as observed experimentally [ 3 ]. In order to check for the finite size effects, we performed our simulations at two different sizes for a typical set of parameters. Fig. 5 shows the rms displacement of small particles versus time on a log-log scale for both low and high free area fractions at the concentrations Pl = 0.30,p2 = 0.49,p3 = 0.21 and Pl = 0.30,p2 = 0.28,p3 = 0.42, respectively. The size of the lattices were 128 x 128 and 256 x 256 with the large particles of size 82 . We see that the data for the different size lattices lie nearly on the same curve and therefore we conclude that the qualitative behavior remains unchanged, due to size of the lattice. Finally, we have also explored the effects of the size of the large molecule (or the ratio of the size of large to small molecules) on the diffusive behavior. For the two sets of free area fractions, at a fixed concentration of the constituents, we have studied the variation of the rms displacement of the small and large particles by varying the size of the large particles. Fig. 6 shows the variation of the rms displacement with time on a l o g - l o g scale. The rms displacement of the small particles does not seem to be affected by the size of the large particles. On the other hand, the variation in the size affects the diffusion of the large particles; the larger the size, the smaller/slower is their diffusion. The largest particles of size 322 has barely begun moving toward the end of our run
R.B. Pandey/Physica A 223 (1996)309-320
318
100
Q
c
E
10
• • 00
8
0
0
~
CO rr
• ••mm~
0.1
lit 1 04
1 0s Time
Fig. 6. RMS displacement of small (open symbols) and large (filled symbols) particles on a 256 x 256 lattice with Pl = 0.30,p2 = 0.28, Ae = 0.60. Large particle of size 82 (circle), 162 (square) and 322 (triangle) are considered. 5-15 independent samples are used.
time. Thus, even in the annealed systems, the large molecules may act as a quenched barriers, at least on a moderate time scale (i.e. t ~ 8 x 104 for our largest particles at high free area fraction Ae = 0.60). A plot of the diffusion constant versus inverse linear size of the large particle is presented in Fig. 7. Note that the diffusion constant does not decrease linearly with the length of the large molecule. In summary, we have presented a simple computer simulation model to study the lateral diffusion in a binary mixture of small phospholipids and large macromolecules. We consider a large mass/size ratio ( > 5 0 ) of the macromolecule and the small particle. We find that the diffusion constant o f the small particle decreases on reducing their concentration at both low and high free area fractions in our annealed system in which both type of particles are mobile. The variation in the diffusion constant in our annealed system seems to be different from that in a corresponding quenched system with large macromolecules fixed. The diffusion constant of the large particles decreases on increasing their concentration, an effect opposite to that observed for the small particles. This finding is consistent with the recent experimental observations on the lateral diffusion in a binary lipid mixture. The decay of the diffusion constant for the small particles on decreasing their concentration seems to suggest the onset of a transition (from slow to fast diffusive regime) with respect to their concentration. The corresponding variation of the diffusion constant of the large particle on increasing their concentration seems to
319
R.B. Pandey/Physica A 223 (1996) 309-320
1.0-
Et~ 0 . 8 rl
--I
"6
0.6-
t6
g o •~ 0 . 4 a
~ 0.2-
o Z
©
0.0
3
0 I
0.04
I
I
I
0.06 0,08 0.10 1/(Linear Size of the Large Particle)
!
0.12
Fig. 7. Normalized diffusion constant of the large particles versus their size with the same statistics as in Fig. 6. show a power law decay. The diffusion constant of the large particles depends on their linear size, and shows a nonlinear decay with their length.
Acknowledgements This work was suggested by H y u k Yu when RBP was visiting his research group at the University o f Wisconsin, Madison. The warm hospitality by H y u k Yu and his research group, in particular Jingwen Ma, is greatfully acknowledged. The author thanks D. Stauffer, G. Ilgenfritz, and Lutz Schlicht for useful discussion and Brad Spencer and Eric Patterson for the help with the local computing set up at the UW-Madison. Some work on the manuscript was done at the Institute for Theoretical Physics ( C o l o g n e University) where the author was visiting as part of his sabbatical. This work is partially supported by a N S F - E P S C o R grant for RBP. The generous support for the computing resources at the Mississippi Center for Supercomputing Research is acknowledged. References [ 1] H.T. Tien, Bilayer Lipid Membrane - Theory and Practices (Marcel Dekker, New York, 1974). [2] T.D. Osborn and P. Yager, Langmuir 11 (1995) 8; A. Tronin, T. Dubrovsky and C. Nicolini, Langmuir 11 (1995) 385; R. Reiter, H. Motschrnann and W. Knoll, Langmuir 9 (1993) 2430.
320
R.B. Pandey/Physica A 223 (1996) 309-320
[3] Hyuk Yu, Private Communication and recent work from his research group, i.e. Langmuir 9 (1993) 1545. [4] K. Tamada, S. Kim and H. Yu, Langmuir 9 (1993) 1545. [5] M.J. Saxton, Biophys. J. 39 (1980) 165. [6] G. Gompper and M. Schick, in: Phase Transition and Critical Phenomena, Vol. 76, C. Domb and J.L. Lebowitz, eds,~(Academic Press, New York, 1994). [7] D. Stauffer, N. Jan, Y. He, R.B. Pandey, D.G. Marangoni and T. Smith-Palmer, J. Chem. Phys. 100 (1994) 6934. [8] J.L. Becklehimer and R.B. Pandey, J. Stat. Phys. 75 (1994) 765. [9] J. Feder, J. Theor. Biol. 87 (1980) 237; R.H. Swendsen, Phys. Rev. A 24 (1981) 504; R.D. Vigil and R.M. Ziff, J. Chem Phys. 93 (1990) 8270; N.M. Svrakic and M. Henkel, J. Phys. (France) I 1 (1991) 791. [ 10] D. Stauffer and A. Aharony, Introduction to Percolation Theory (Taylor and Francis, London, 1992); M. Sahimi, Application of Percolation (Taylor and Francis, London, 1994). [11] A. Drory, I. Balberg and U. Alon, Phys. Rev. A 43 (1991) 6604. [ 12] K. Binder, ed., Application of Monte Carlo Simulation in Statistical Physics (Springer, Berlin, 1986). [13] S. Havlin and D. Ben-Avraham, Adv. Phys. 36 (1987) 695. [ 14] R.B. Pandey, D. Stauffer, A. Margolina and J.G. Zabolitzky, J. Stat. Phys. 34 (1984) 427. [15] R.B. Pandey, Physica A 187 (1992) 77. [16] R. Kutner and K. Kehr, Phil. Mag. A 48 (1983) 188. [ 17] G. Ilgenfritz and E Runge, Physica A 181 (1992) 69. [18] M.E Allen and D.J. Tildesley, Computer Simulation of Liquids (Clarendon, Oxford, 1989).
Physica A 223 (1996) 309-320
Lateral diffusion in a binary lipid system by a computer simulation model R.B.
Pandey
Department of Physics and Astronomy, The Program in Scientific Computing, University of Southern Mississippi, Hattiesburg, MS 39406-5046, USA
Received 4 October 1995; revised 19 October 1995
Abstract A computer simulation model is used to study the lateral diffusion in a two-dimensional binary mixture of large and small particles. The small particle represents the phospholipid while the large particle (square) models the planar macromolecule in the binary mixture. The variation of the rms displacement with time of these molecular species are studied as a function of their concentration and their diffusion constants are examined. The concentration of the large particles (pl), small particles (p2), and the vacancy (P3) are varied at the relative area fraction Ae = p3/(P2 q- P3) = 0.30 (low) and 0.60 (high) free area fractions. The diffusion constant of the small particles decays on decreasing their concentration - this decay signals the possibility of a transition between a fast diffusive regime (at high concentration) and a slow diffusive regime (at low concentration p2). The diffusion constant of the large molecules on the other hand shows an opposite behavior, it decays linearly on increasing their concentration (pl). The diffusion constant of the large molecule decays with their size while that of the small particles remains less affected by varying the size of the large molecules. The decay of the diffusion constant of the large molecule with its linear size is found to be nonlinear.
1. I n t r o d u c t i o n The molecular mobility of phospholipids in cell membranes plays a very important role in governing their physical and chemical properties [ 1 ]. Of our particular interest is the lateral diffusion of lipids in biphasic phospholipid monolayers [2,3]. The system consists of a binary mixture of two species - small phospholipids and large lipid macromolecules with their weight ratios of the order of 50:1 to 100:1. Extensive experimental studies [3,4] such as surface pressure measurements with FRAP method and fluorescence microscopy have been recently carried out to measure the diffusion constant of such binary mixtures as DLPC (L-alpha-Dilaroylphosphatidylcholine). Several 0378-4371/96/$15.00 (~) 1996 Elsevier Science B.V. All rights reserved SSDI 0378-4371 (95)00366-5
310
R.B. Pandey/Physica A 223 (1996) 309-320
interesting findings have been reported on the variation of the diffusion constant with the concentration of the constituents, their size, and the nature of the interfacial substrate among other parameters. The mobility of the two species obviously depends on several factors such as concentration of these molecules, free surface/volume, their relative size, shape, interaction, and temperature, etc. Understanding the diffusion in such a complex mixture is one of the most difficult theoretical problems due to strong correlations from small to large scales. Even the approximate theoretical methods such as mean field theories including selfconsistent field and effective medium [ 5] approaches are limited to monodisperse systems; most of these methods are known to be inadequate for taking into account the effects of correlations. Therefore, computer simulation methods remain one of the very important tool to understand the diffusion in such complex systems. Computer simulation models have been recently used to study the phase separation and phase transition in simplified amphiphilic systems [6,7]. The diffusion of amphiphiles are particularly used to monitor the evolution of micelle formation as the system evolves through a bicontinuous phase in a model of micoemulsion [7]. Because of the limitations on the computing resources (memory and speed), it is not feasible to take into account all the relevant details of such a binary system. Therefore, we will limit ourselves to simplified model. We study here the diffusion in a binary lipid system consisting of planar macromolecules and a small phospholipid by a simplified computer simulation model.
2. Model
We consider a discrete L x L square lattice. The small phospholipid molecule is modelled by a point-like particle while the large lipid macromolecule by a solid square of size l x I. Thus, a small particle occupies a lattice site and a large particle 12 lattice sites. We would like to initialize the lattice by randomly distributing these particles with a concentration Pl of large, and P2 of small particles; the concentration P3 of vacancy (i.e. empty sites) is P3 = 1 - p l - P 2 . A technical problem regarding the random packing arises when we attempt to distribute the large particles beyond a certain concentration (Plc), known as the first jamming/saturated coverage [8,9] or percolation threshold [ 10]; generally Plc depends on the size and shape of the objects, and decreases with the size of the object with power-laws for several well-defined shapes [ 8,11 ]. However, in order to get around this problem, one may distribute the large particles in a regular fashion, and then move them around randomly for a predetermined period of time before quenching them to obtain their random distribution. Such process of annealing followed by quenching may not provide a complete random distribution, but offers an alternate solution to jamming coverage. Our system is close to an annealed system (see below), therefore, we believe that this method of initializing our lattice system does not affect our results. Thus, we first distribute large particles on a fraction Pt of the lattice sites and then a fraction P2 of the lattice sites are randomly occupied by small particles. An excluded volume (hardcore) interaction is considered between the particles; a large
R.B. Pandey /Physica A 223 (1996)309-320
311
particle occupies l 2 lattice sites while a small particle occupies a lattice site. A lattice site cannot be occupied by more than one particle at a time. We stir the system by moving these particles randomly for a period of time as follows. For small particles, we choose a particle and one of its neighboring site randomly; if the neighboring site is empty then we move the particle. However, if the neighboring site is occupied by a small particle or by part of a large particle, then the particle remains at its old position. Essentially the same procedure is applied in moving a large particle as well. We attempt to move all the sites of an entire macromolecule by one lattice constant in a randomly chosen direction. Thus, for a successful move, there must be at least l neighboring sites vacant in the direction of the move. The probability of finding l empty neighboring sites along the edge of the large particle in the direction of the hopping is usually small in a dense system. The motion of the large particles is therefore expected to be slow. Furthermore, even though the center of mass of the large particle is shifted by one lattice constant in a successful hopping, l sites simultaneously become vacant and occupied in moving a large particle. Therefore, we attempt to move the large particle with a hopping rate ~ 1 / l (see below). The attempt to move each particle regardless of its success to hop to its neighboring site/sites is considered in defining the time step. The definition of Monte Carlo Step is generally arbitrary [ 12] ; an attempt to move each particle once is usually defined as one Monte Carlo Step (mcs). In our case here, there are two particles of different sizes, therefore we define a unit time step in which each small particle attempts to move once while a large particle moves with a probability 1 / l . A periodic boundary condition is used to move the particles across the boundaries. We consider two systems. ( 1 ) A quenched system, in which we move all the particles in the system for a pre-assigned period of annealing time, and then we quench the large particles. The "so-called" random distribution of large particles thus achieved in such an annealing process acts as a quenched random barrier for the mobility of the small particles. Such quenched porous media may be compared to a percolating porous media. However, one should be careful in comparing the physical properties of such porous medium with that of ordinary site/bond percolation, since not only the percolation threshold of such systems are different but also the morphology, and, therefore, the geometrical correlations may differ substantially. (2) An annealed system, in which we let the particles move around for a sufficiently long time so that we achieve their asymptotic transport behavior. The whole simulation is repeated for a number (N~un) of independent samples to obtain a reliable average of the transport quantities. During the simulation we keep track of the displacement of all particles and their center of mass, and we evaluate the rms displacement of each particle periodically during the entire course of time. The leading power-law behavior of the rms displacement (R) of a particle with time t is described by R = 2Dt k ,
(1)
where k is a power-law exponent, and 2D is a prefactor. For most annealed systems, k = 1/2 in the asymptotic regime, and the stochastic motion of the particles is called
R.B. PandeylPhysica A 223 (1996)309-320
312
,,,,/
(a)
10098" 76-
E
s-
~ 131
4-
~
3-
8
10-
| I
I
i
1 04
1 0s Time
I
(b) 2.22 -
2.20 -
aI
c
°
2.18-
A All ~
ffJ
& j
2.15-
All
n-
•
Vo°
•.
V..,~ .V~ V ...V
AI !re ~1 "~o
2.14-
I
* •
•
Vo
V o
~o 2.10-
6x 10 4
1 05 Time
F i g . 1.
R.B. Pandey/Physica A 223 (1996) 309-320
313
Fickian or diffusive with the Einstein diffusion constant D. The motion o f a particle can, however, be anomalous and nondiffusive in many systems such as in random percolating system at percolation threshold [10,13,14] or in a complex mixture o f interacting systems [ 15] where the exponent k v~ 1/2. In the latter case, the "diffusion constant" terminology for the prefactor ( D ) may be misleading. We will, however, stick to the diffusion constant [16,17] since the motion is diffusive ( k = 1 / 2 ) for most o f the concentration regime considered here. The diffusion constant D depends on the heterogeneity o f the host space [ 14] and the concentration o f the mobile species [ 16] which is the main focus o f our study in the following. We w o u l d like to mention that there are different methods o f estimating the diffusion constant. In continuum simulations, particularly in Molecular Dynamics, one usually evaluates the velocity autocorrelation function, and averaging over time then gives an estimate o f diffusion constant [ 18]. Alternatively, since the motion is diffusive, the instantaneous derivative o f the average value o f R 2 with time will give us the diffusion constant which should approach a constant value in the asymptotic limit ( t ~ cx~). The intercept o f a l o g - l o g plot o f the rms displacement versus time may also give us an estimate o f the diffusion constant which may vary somewhat with the number o f data points particularly when the data are fluctuating. In order to avoid confusion, we will concentrate to an effective diffusion constant ( 2 D ) , which is obtained by evaluating R 2 / t , where t is the largest time in our simulations. Since we are interested mainly in the variation o f the diffusion constant with the concentration o f the constituents rather than its absolute value, the evaluation methods should not affect our findings. Therefore, we will use the simple method here. Furthermore in order to reduce the relative fluctuations and to maintain the consistency, the period o f simulation and time t at which the diffusion constant is evaluated is kept the same for all the simulations to study its relative variation with the concentrations o f the two types o f particles.
3. Results and discussion The simulations are performed on square lattices o f size 128 x 128 and 256 × 256 with various sizes o f the large particles, 82 to 322 . The large particles are used to model the effects o f large macromolecules and their m a s s / s i z e reflects a relative measure with respect to small particles which represent the small phospholipids. The number o f independent runs Nrun are limited and varied with Nrun = 5 - 2 0 in order to optimize the computing resources and to reduce the statistical fluctuations. As we mentioned
Fig. 1. RMS displacement of small particles versus time on a log-log scale with (a) on the full time scale while (b) is only in the long time regime. Lattice of size 128 x 128 was used with the large particles of size 82 with 20 independent samples. The relative free area fraction Ae 0.30 at the concentrations Pl = 0.1,p2 = 0.63 (circle), Pl = 0.2,p2 = 0.56 (square), Pl = 0.3,p2 = 0.49 (triangle), Pl = 0.4,p2 = 0.42 (inverted triangle), Pl = 0.5,p2 = 0.35 (parallelogram). Open symbols represent the rms displacement of the particles in quenched system while the filled symbols represent the rms displacement in the annealed system. =
R.B. Pandey/Physica A 223 (1996) 309-320
314
0
1.00
0.98 o 1= n
0
0.96 -
0,94 -
8 0.92 a
) )
0.90
o Z
O 0~88 •
0.0
I 0.1
I 0.2 Concentration
I 0.3
• 0
of Small
0.4
I 0.5
I 0.6
Particles
Fig. 2. Normalized diffusion constant of small particle versus their concentration (P2) at the free area fraction Ae = 0.30. The statistics is the same as in Fig. 1. The open circle represent the diffusion constant of the particle in the quenched system while the filled circle, in the annealed system. The diffusion constants are normalized by their maximum value.
before, we consider two systems for comparison. (1) The quenched system in which the large macromolecules are quenched and act like quenched barriers for the mobile phospholipids. (2) The annealed system, in which both large and small molecules are mobile with hopping attempt rates inversely proportional to their linear size. We limit ourselves to two values of the free area fraction, Ae = P2/(Pz + P3), a low free area (Ae = 0.30) and a relatively high free area (Ae = 0.60). For both sets of free area fractions, the concentration of the large macromolecule (Pl), small phospholipids (P2), and the free area (P3) are varied with Pl + P 2 q-P3 = 1. All the data are generated on a Cray YMP and a Cray J916 machines; about 100 hours of CPU are used. We have studied the variation of the displacement and the root mean square (rms) displacements of both small and large particles and their center of mass. We monitor the displacements primarily as an alternate method to check whether the system has reached the asymptotic regime. The variation of the rms displacement with time is the main quantity we analyse here. The large particles move relatively slow, therefore we will mainly concentrate on the rms displacement of the small particles except toward the end to discuss the effects of the size of the large particles on their transport. The variation of the rms displacement of the small particle with time, at a free area fraction, Ae = 0 . 3 0 , is presented in Fig. 1 for both quenched and annealed systems for various concentrations of the constituents. From a first glance on Fig. la, it appears that there is no effect on
R.B. Pandey/Physica A 223 (1996) 309-320 1.00
315
-
•
0
©
g
© 0.95
-
II.
0
t~
E
0.90-
"6 c
0.850
i5 0 . 8 0 lo
g
o 0.75 z
0.70
0.00
I
I
I
0.05
0.10
0.15
Concentration
• I
I
I
I
0.20
0.25
0.30
0.35
of S m a l l
Particles
Fig. 3. Normalized diffusion constant versus concentration of the small particles, similar to Fig. 2, for the free area fraction Ae = 0.60. The statistics is the same as in Fig. 2.
changing the concentration ( p 2 ) of the small particles with a little difference in the temporal variations of the rms displacement for the quenched and the annealed systems. All the data seem to collapse nearly to a single curve with a linear diffusive power-law fit, i.e. R N t l / 2 especially in the long time (asymptotic) regime; from Eq. (1), k = 1/2. Even the intercept of these data points look nearly the same in Fig. la. However, a close examination of these curves (Fig. lb) reveal the fact that the diffusion constant (i.e. the intercept) depends on the concentration of the constituents; diffusion constants at some concentrations in the annealed system differ from those in the quenched systems. As we mentioned before, there are several methods to evaluate the diffusion constant. We estimate an effective diffusion constant De = 2D = ( R 2 ) / t as t ---+~ , the largest time in our study. In Fig. 2 we present the variation of the effective diffusion constant of the small particles with their concentration at the low free area fraction (Ae = 0.30). As expected, the diffusion constant for the small particles in the annealed system has comparatively less fluctuations than the corresponding values in the quenched system. In the annealed system, the diffusion constant decays on reducing the concentration of the small particles until a characteristic concentration (P2c) below which it remain nearly a constant to a low value. This may be a signal of a transition, from a relatively fast diffusion at/)2 > P2c to a slow diffusion at p2 < p2c. This transition may be different from the conductor-insulator transition in a static percolation [7]. The former is a dynamic heterogeneous system while the latter is a quenched inhomogeneous
R.B. Pandey/Physica A 223 (1996) 309-320
316 t.0-
0h91:: e5 ¢1 ¢D
O
0.8o c e~
0.7-
8
g 0.6
ca N
0.5
O
z
0.4
I
0.2
I
I
I
0.4 0.6 0.8 Concentration of Large Particles
I
1,0
Fig. 4. Normalized diffusion constant versus concentration of the large particles at the free area fraction Ae = 0.30 (open circle) and Ae = 0.60 (filled circle). The statistics is the same as the previous figures.
system. Unfortunately, a considerable amount of more computing resources are needed to quantify the possibility o f such a transition in our annealed system. The diffusion constant in our quenched system, on the other hand, shows a peak before falling off at nearly the Same value o f P2 as the diffusion constant i n the annealed system. However, since there is m o r e fluctuations in the quenched case, we should be Cautious in such conclusions for the quenched system. We have carried out a similar investigation also at a relatively large free area fraction (Ae = 0 . 6 0 ) . The variation o f the corresponding diffusion constant of the small particles with their concentration is shown in Fig. 3. It is relatively harder to evaluate the diffusion constant more accurately due to large fluctuations in the data here than that with low free area fraction. However, as with the low free area fraction, the data with the annealed system shows less fluctuations than that in the quenched system. We see a decay in the diffusion constant on reducing the concentration (P2); the peak at the concentration p2 "~ 0.30 may be an artifact o f fluctuations. The trend in the variation of the diffusion constant with the concentration in the quenched system is even less predictable which, nevertheless, seems to suggest an opposite trend to that in the annealed system. The diffusion o f large molecules is studied in the annealed system. The variation o f the diffusion constant with their concentration is presented in Fig. 4 at the free area fractions Ae = 0.30 and 0.60. Except at the low concentrations (pl = 0.1 - 0.4) at Ae = 0.30, the diffusion constant seems to decay linearly on increasing the concentration
R.B. Pandey/Physica A 223 (1996) 309-320
317
n
100E
9-
"6
B-
"E
~-
OU
7-
E
8
6-
== co
O
n:
[]
•
Q
0 1 04
1 05 Time
Fig. 5. RMS displacementof small particle versus time on a log-log scale at the free area fraction Ae = 0.30 (circle) and Ae = 0 . 6 0 ( s q u a r e ) . The size of the large particle is 82 in all the runs. Lattice of size 128 x 128 (open symbols) is used with 20 independentsamples while the lattice of size 256 x 256 (filled symbols) is used with 5 and 15 independentruns for Ae = 0.30 and 0.60, respectively. of large particles. Note that the diffusion constant (D1) of the large particle shows an opposite trend to that for the small particle as observed experimentally [ 3 ]. In order to check for the finite size effects, we performed our simulations at two different sizes for a typical set of parameters. Fig. 5 shows the rms displacement of small particles versus time on a log-log scale for both low and high free area fractions at the concentrations Pl = 0.30,p2 = 0.49,p3 = 0.21 and Pl = 0.30,p2 = 0.28,p3 = 0.42, respectively. The size of the lattices were 128 x 128 and 256 x 256 with the large particles of size 82 . We see that the data for the different size lattices lie nearly on the same curve and therefore we conclude that the qualitative behavior remains unchanged, due to size of the lattice. Finally, we have also explored the effects of the size of the large molecule (or the ratio of the size of large to small molecules) on the diffusive behavior. For the two sets of free area fractions, at a fixed concentration of the constituents, we have studied the variation of the rms displacement of the small and large particles by varying the size of the large particles. Fig. 6 shows the variation of the rms displacement with time on a l o g - l o g scale. The rms displacement of the small particles does not seem to be affected by the size of the large particles. On the other hand, the variation in the size affects the diffusion of the large particles; the larger the size, the smaller/slower is their diffusion. The largest particles of size 322 has barely begun moving toward the end of our run
R.B. Pandey/Physica A 223 (1996)309-320
318
100
Q
c
E
10
• • 00
8
0
0
~
CO rr
• ••mm~
0.1
lit 1 04
1 0s Time
Fig. 6. RMS displacement of small (open symbols) and large (filled symbols) particles on a 256 x 256 lattice with Pl = 0.30,p2 = 0.28, Ae = 0.60. Large particle of size 82 (circle), 162 (square) and 322 (triangle) are considered. 5-15 independent samples are used.
time. Thus, even in the annealed systems, the large molecules may act as a quenched barriers, at least on a moderate time scale (i.e. t ~ 8 x 104 for our largest particles at high free area fraction Ae = 0.60). A plot of the diffusion constant versus inverse linear size of the large particle is presented in Fig. 7. Note that the diffusion constant does not decrease linearly with the length of the large molecule. In summary, we have presented a simple computer simulation model to study the lateral diffusion in a binary mixture of small phospholipids and large macromolecules. We consider a large mass/size ratio ( > 5 0 ) of the macromolecule and the small particle. We find that the diffusion constant o f the small particle decreases on reducing their concentration at both low and high free area fractions in our annealed system in which both type of particles are mobile. The variation in the diffusion constant in our annealed system seems to be different from that in a corresponding quenched system with large macromolecules fixed. The diffusion constant of the large particles decreases on increasing their concentration, an effect opposite to that observed for the small particles. This finding is consistent with the recent experimental observations on the lateral diffusion in a binary lipid mixture. The decay of the diffusion constant for the small particles on decreasing their concentration seems to suggest the onset of a transition (from slow to fast diffusive regime) with respect to their concentration. The corresponding variation of the diffusion constant of the large particle on increasing their concentration seems to
319
R.B. Pandey/Physica A 223 (1996) 309-320
1.0-
Et~ 0 . 8 rl
--I
"6
0.6-
t6
g o •~ 0 . 4 a
~ 0.2-
o Z
©
0.0
3
0 I
0.04
I
I
I
0.06 0,08 0.10 1/(Linear Size of the Large Particle)
!
0.12
Fig. 7. Normalized diffusion constant of the large particles versus their size with the same statistics as in Fig. 6. show a power law decay. The diffusion constant of the large particles depends on their linear size, and shows a nonlinear decay with their length.
Acknowledgements This work was suggested by H y u k Yu when RBP was visiting his research group at the University o f Wisconsin, Madison. The warm hospitality by H y u k Yu and his research group, in particular Jingwen Ma, is greatfully acknowledged. The author thanks D. Stauffer, G. Ilgenfritz, and Lutz Schlicht for useful discussion and Brad Spencer and Eric Patterson for the help with the local computing set up at the UW-Madison. Some work on the manuscript was done at the Institute for Theoretical Physics ( C o l o g n e University) where the author was visiting as part of his sabbatical. This work is partially supported by a N S F - E P S C o R grant for RBP. The generous support for the computing resources at the Mississippi Center for Supercomputing Research is acknowledged. References [ 1] H.T. Tien, Bilayer Lipid Membrane - Theory and Practices (Marcel Dekker, New York, 1974). [2] T.D. Osborn and P. Yager, Langmuir 11 (1995) 8; A. Tronin, T. Dubrovsky and C. Nicolini, Langmuir 11 (1995) 385; R. Reiter, H. Motschrnann and W. Knoll, Langmuir 9 (1993) 2430.
320
R.B. Pandey/Physica A 223 (1996) 309-320
[3] Hyuk Yu, Private Communication and recent work from his research group, i.e. Langmuir 9 (1993) 1545. [4] K. Tamada, S. Kim and H. Yu, Langmuir 9 (1993) 1545. [5] M.J. Saxton, Biophys. J. 39 (1980) 165. [6] G. Gompper and M. Schick, in: Phase Transition and Critical Phenomena, Vol. 76, C. Domb and J.L. Lebowitz, eds,~(Academic Press, New York, 1994). [7] D. Stauffer, N. Jan, Y. He, R.B. Pandey, D.G. Marangoni and T. Smith-Palmer, J. Chem. Phys. 100 (1994) 6934. [8] J.L. Becklehimer and R.B. Pandey, J. Stat. Phys. 75 (1994) 765. [9] J. Feder, J. Theor. Biol. 87 (1980) 237; R.H. Swendsen, Phys. Rev. A 24 (1981) 504; R.D. Vigil and R.M. Ziff, J. Chem Phys. 93 (1990) 8270; N.M. Svrakic and M. Henkel, J. Phys. (France) I 1 (1991) 791. [ 10] D. Stauffer and A. Aharony, Introduction to Percolation Theory (Taylor and Francis, London, 1992); M. Sahimi, Application of Percolation (Taylor and Francis, London, 1994). [11] A. Drory, I. Balberg and U. Alon, Phys. Rev. A 43 (1991) 6604. [ 12] K. Binder, ed., Application of Monte Carlo Simulation in Statistical Physics (Springer, Berlin, 1986). [13] S. Havlin and D. Ben-Avraham, Adv. Phys. 36 (1987) 695. [ 14] R.B. Pandey, D. Stauffer, A. Margolina and J.G. Zabolitzky, J. Stat. Phys. 34 (1984) 427. [15] R.B. Pandey, Physica A 187 (1992) 77. [16] R. Kutner and K. Kehr, Phil. Mag. A 48 (1983) 188. [ 17] G. Ilgenfritz and E Runge, Physica A 181 (1992) 69. [18] M.E Allen and D.J. Tildesley, Computer Simulation of Liquids (Clarendon, Oxford, 1989).