摘要:

Simulating the self-assembly of model membranes Maddalena Venturoli and Berend Smit Department of Chemical Engineering, Universiteit van Amsterdam, Nieuwe Achtergracht 166, 1018 WV Amsterdam, The Netherlands Received 9th August 1999, Accepted 9th September 1999, Published 15th September 1999 Dissipative particle dynamics simulations are presented of the self assembly of surfactant bilayers. The effect of changes in the chain length and stiffness of the surfactants on the properties of the model membranes are studied. We observe that changes of the stiffness have significant effects if these changes are made close to the head group of the surfactant. If, on the other hand, changes are made at the end of the tail of the surfactant, the properties of the bilayer are similar to the properties of a bilayer consisting of flexible chains. 1 Introduction The role of the lateral pressure profile in biological membranes has been the topic of discussions in the literature in the context of the structure and function of membrane proteins1 or in the mechanism of anesthesia.2 An important aspect in these discussions is that there are no experimental methods to determine a pressure profile.At present, we have to rely on simulations or theory.2 By using molecular simulation it is, in principle, possible to study an all atom model of a membrane (see for example ref. 3 and references therein). However, to do so one has to start already from a reliable initial configuration of the membranes, else it would require far too much cpu time to observe the self-assembly of such a model membrane. One can also adopt a coarse-grained approach in which a less detailed model of a lipid bilayer is used.Goetz and Lipowsky4 have shown that with such a model one can study the self-assembly of a (model) membrane using conventional molecular dynamics. In this work, we introduce an alternative method to simulate the self-assembly of a model membrane, which is based on the dissipative particle dynamics (DPD) simulation technique.5,6 Furthermore, we use a potential function that is much ’softer’ than the Lennard-Jones potential conventionally used in MD simulations. By making these choices we are able to use a time step that is significantly larger (by a factor of 2 to 5) than the one used in conventional molecular dynamics simulations.Goetz and Lipowsky have also shown that the use of periodic boundary conditions and a fixed number of surfactants per area leads to a bilayer that has a given value of the surface tension. Biological membranes, however, have a state which is essentially tensionless. To ensure such a tensionless state Goetz and Lipowsky performed several simulations to determine the area per bilayer lipid that gives a state of zero tension. Here, we show that we can mimic the experimental situation by a constant surface tension simulation.7 Similar to constant pressure simulation, we impose a zero value for the surface tension and from the simulation we obtain the average area per surfactant.Such simulation has the advantage that we do not have to perform several—relatively expensive—simulations to locate the area of zero tension. PhysChemComm, 1999, 10 As an illustration of the type of simulations that can be performed with our method, we present some results in which we study the effect of the chain length and stiffness of our model lipids on the properties of the bilayer. c 2 Simulation techniques In DPD the forces are composed of (soft) repulsion conservative forces, pairwise dissipation forces, and pairwise random forces. The force acting on a particle i is then given by: , (1) F + F + F ) R ij D ij ( ijC i å�i j D , r� )( (r ×v )r�ij ij ij ij r� fFF ij ij r r - r r� ) (rij ij c ij = with D ijR ijC ij ( rij < ) rc 2 2 - = hwR rij ) (R r , )] ( 2 = sw where rij = rj – ri is the distance vector between particles i and j, rij =|rij| and vij = vj – vi.For our simulations we have used the soft-core repulsive force that has been used in many other DPD simulations:8 ì (1 a0 ï ij F = íïî According to Español and Warren6 the DPD technique samples a Boltzmann distribution with a potential related to the conservative force: Fw -ÑU provided that the weights wD and wR and the amplitudes s and h of the dissipative and random forces [eqn. (2)] satisfy the following relations: w r) (= C D = [ � ) BT k s = h where wR(r) is taken as in ref. 8. As shown by Willemsen et al.,9 the existence of a Hamiltonian implies that DPD can be combined with (2) (3) (4)Monte Carlo methods. In this work, we combine DPD with a Monte Carlo scheme to ensure that the simulations are performed at constant surface tension.In an ensemble of constant temperature T (or constant b = 1/kBT), number of particles, N, and surface tension, g, the area of the bilayer, A, is allowed to change in such a way that the total volume, V, of the box is constant. The partition function of such an ensemble can be written as: 1 N N A Q r ) U( ( d - b = l dr exp 3N N! L 2 3N 0 = 3N N! 3N 0 - g [ ] 2N N N U (5) d (s ) s exp 2N ò s exp lld b - - b [ ] N U(s ( ò ò d exp l L3NN! z ò l d exp y x 0 N bgl [ ] N 2 ò in which we have introduced the scaled coordinates, s, and the order parameter, l , by defining ) ) , l µ - l® l = ¢ bgl [ ] [ ]) ) N 2 s ] (6) l2 where L0 is an arbitrary unit of length.This form of the order parameter ensures that the total volume remains constant. The probability of finding a configuration with positions sN and l, is given by: ( [U - b gl { } To simulate at a constant surface tension we perform, during the simulation, Monte Carlo moves in which we change the order parameter l. From eqn. (6) the probability of accepting such a move is given by ¢N b[ ( U N gl¢2 2 s ) s ) exp{- exp{ ]} ]} U( [ -- b - gl To ensure microscopic reversibility, our simulations are done in cycles.In each cycle we perform a random number of time steps of ordinary DPD simulations and an attempt to change l . On average we perform 50 time steps (Dt = 0.06) per cycle. To test this algorithm we have compared several simulations with different initial surface areas per surfactant molecules. Within 100–300 Monte Carlo cycles an equilibrium surface area is reached independent of the initial area. To ensure that the bilayer is always formed in the x,y plane, we use a rectangular box in which the x,y plane has a significantly smaller area than the x,z and z,y planes. Given the constraint on the number of surfactant particles a bilayer can only form in the x,y plane. The formation of the bilayer is done using conventional DPD. Once the bilayer is stable, we continue the equilibration with the combined DPD and zero-surface tension scheme till a constant area per head group was achieved.This equilibration took about 5000–10000 cycles depending on the type of surfactant. After equilibration the production runs were at least 20000 cycles. =L = rN s ( acc( ò ls , s s l , exp ) 3 Model In our simulations, we distinguish three types of particles that model water (w) and the hydrophilic head (h) and hydrophobic tail (t) segments of the surfactants. All interactions are described with the conservative force given by eqn. (4) with parameters chosen in order to mimic the hydrophobic and hydrophilic interactions of the bilayer lipids with water: (7) aww = ahh = att = 25 aht = 40 awt = 30 We have used reduced units. The unit of mass is defined by specifying that the mass of all DPD particles is 1, the length scale is defined by the cut-off of the potential (rc = 1), and the temperature scale is defined by the potential in such a way that a reduced temperature of 1 corresponds to T= a/kB, where a is the unit of energy in which the parameters of eqn.(7) are expressed. In addition, surfactant molecules are constructed by connecting h and t atoms via harmonic springs: (8) 1 i i, Uspring (r + +1 i i, r r k 2 2 ) 2 ( ) Ubend 0 i i 1, i r ( - , r + 2 × r k i+1,i q ö÷ = with kr = 2. This potential, together with the repulsive force of eqn.an average bond-length of 0.63.A surfactant molecule consists of a linear chain of h and t units. For example, a linear chain with one head group and seven tail units is denoted by ht7. The effect of chain length on the bilayer's properties can be studied by changing the number of tail units. Unsaturated carbon bonds can change the stiffness of a membrane lipid. We model the stiffness of the chain by introducing a bond bending potential: i 1, k q 2 ) - 0 ÷ (9) 2 i 1, i i 1, r + - æ ri-1,i i ççè r - where q is the angle between two consecutive bonds. We have used k q = 3 and q0 = p . In the notation we indicate the presence of the bond bending potential by a capital T for the atom i in eqn. (9), for example, ht2Tt indicates a chain of 5 beads in which the last three atoms of the tail have the bond bending potential. It is important to note that both the spring and bond-bending potential are part of the conservative forces for the DPD program.All simulations have been performed at k surfactants. The concentration, c defined as ø BT = 0.2, with 100 s, of the surfactants, cs == + Nt + N + Nt Nh Nw h = was equal to approximately 0.39 for all the simulations. Where Nh, Nt and Nw are, respectively, the number of head, tail and water particles. The volume of the box was chosen to have a density of about 2.9 (at such density and concentration of surfactants where a bilayer was formed). 4 Results and discussion The self-assembly of the bilayer is shown in the movie of Fig.1. We start with a random distribution of surfactants(see Fig. 2). In the first stage these surfactants form a cluster. This cluster then takes the shape of a cylinder-like structure. This cylindrical micelle is not stable and shows large fluctuations in shape. Some of these fluctuations result in a ’percolation’ of the micelle across the periodic boundary conditions and a bilayer is formed (see Fig. 2). Although it is interesting to observe the self-assembly of a bilayer, the mechanism will be very much dependent on the system size. Fig. 1 Animation of the self-assembly of a bilayer starting from a random distribution of surfactants. A snapshot is taken every 200 time steps and the movie is composed of 170 snapshots. The total number of surfactants is 100 and there are 900 water particles.The water particles are not shown in the movie. The initial frame is a projection on the side plane x,z (where z is the horizontal axis), and during the movie the system rotates around the 'up–down' axis x. Fig. 2 The initial (left figure) and final (right figure) configurations of the simulation of the self-assembly. The total number of surfactants is 100 and there are 900 water particles. The water particles are not shown. The bilayer is in the x,y plane. The effect of the tail length and stiffness on the average area per surfactant is shown in Fig. 3. The area per surfactant was calculated by computing the average area during the simulation and dividing this by half the number of surfactants. In all our simulations the number of surfactants was equal in both layers during the entire simulation.As the tail length increases, the area per surfactant increases. It is interesting to compare these results with the theoretical calculations of Cantor10 on a lattice model. To make this comparison we used the following mapping of parameters. Since our DPD particles represent a collection of atoms which are more flexible than Cantor’s lattice model, it is reasonable to represent 2 lattice sites by one DPD particle. To ensure similar length scales in both studies, we have scaled the areas to give the same area per surfactant for chains with length 7. Fig. 3 shows that both studies predict that the area increases linearly with chain length. Considering the ad-hoc character of our mapping it would be too preliminary to draw conclusions on the very good agreement between these studies.In Fig. 3 we also compare the effect of the stiffness on the area per surfactant. This figure shows that the effect depends on the location of the bond bending potential. If the bond bending potential is at the end of the tail, the effect on the area is very small and the results are nearly identical to a completely flexible chain. However, if we introduce a stiffness very close to the head group of the surfactant, we find a pronounced decrease of the area per surfactant compared to flexible surfactants. The stiffness close to the head groups gives a stronger ordering of the molecules and as a result of this the mutual interactions result in a more compact structure.Fig. 3 Effect of the surfactant structure on the area per surfactant. In the top figure we compare three types of surfactant completely flexible htn, surfactants with a stiffness at the end of the tail htn–2Tt, and surfactants with a stiffness at the head hTtn–1. In the bottom figure we compare our simulation results with the calculation of Cantor for completely flexible chains (for the definitions of the area and chain length see text). The pressure and density profiles for a flexible surfactant of 6 beads are presented in Fig. 4. We define the pressure profile as the local surface tension: z) ( z) ( g( ) z = t pn - p where pn(z) and pt(z) are the normal and tangential components of the pressure tensor.Español and Warren6 have shown that the ensemble averages can be taken considering only the conservative force. Therefore in the calculation of the pressure tensor we consider the contribution of those forces only. The normal componentdoes not depend on z. Since we simulate at zero surface tension, the integral of this pressure profile is zero. We find close to the head group first a minimum in the pressure profile followed by a large maximum if we move into the interior of the bilayer. In the middle of the bilayer we find a maximum again. The shape of the pressure profiles looks very similar to the ones presented by Goetz and Lipowsky.4 The mean field calculation of Cantor focuses on the pressure profile in the region of the tails, which roughly corresponds to z Î [–1.5,1.5] in Fig.4. Our results in this region are in good (qualitative) agreement with Cantor's. Fig. 4 Density r(z) (top figure) and pressure profile g(z) (bottom figure) as functions of the distance from the middle of the bilayer z for 100 ht6 surfactants. w gives the density profile of the water particles. The effect of changes of the chain length of the surfactants on the pressure profile is shown in Fig. 5. We observe that the maximum close to the head group increases significantly with chain length while the depths of the minima remain nearly constant. However, the minima close to the middle of the bilayer broaden to ensure that the area of the pressure profile remains zero.Fig. 5 Pressure profile g(z) as a function of the distance z of the average position from the head group for 100 htn surfactants (only half of the pressure profile is shown). The effect of the stiffness of the surfactants on the pressure profile is shown in Fig. 6 for short and long chains. As can be expected from the results for the area per surfactant, if we increase the stiffness at the tail of the surfactant there is very little effect, no matter how long is the tail. However, if we increase the stiffness close to the head group of the surfactant we observe a significant effect: the minimum close to the head group gets deeper and the maximum increases compared to a flexible chain. It is interesting to note that this effect is stronger for the shorter chain.A similar effect is observed by Cantor.10 In ref. 10 the effects of unsaturated bonds was studied by imposing a preference for a 90 degree angle between three subsequent sites compared to a 180 degree angle for a saturated chain. Also Cantor found that close to the head group this change in the surfactant conformation has the largest effect. Fig. 6 Pressure profile g(z) as a function of the distance z of the average position from the head group. In the left figure we compare chains of 6 beads in which we varied the location of the bond bending potential and in the right figure chains of 10 beads. Only half of the pressure profile is shown. 5 Concluding remarks In this work, we have shown that DPD type of simulations can be used to simulate the self-assembly of bilayers.We found it convenient to use a zero surface tension ensemble since this avoids the need of having to perform several different simulations to compute the surface area for which the surface tension is zero. We have also presented the results of some preliminary calculations in which we study how changes in the structure of the surfactants affect the physical properties of the bilayer. An interesting conclusion is that stiffness has a large effect on the properties, provided that this stiffness is located close to the head group of the surfactants. Acknowledgements We would like to thank Drs A. Killian, B. de Kruijff, E. van den Brink, M. M. Sperotto and R. Cantor for many stimulating discussions. This work is part of the research program of the 'Stichting voor Fundamentel Onderzoek der Materie (FOM)', which is financially supported by the 'Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)'. Additional support by the Netherlands Organization for Scientific Research (NWO) through PIONIER is acknowledged. References 1 B. de Kruijff, Nature, 1997, 386, 129. 2 R. S. Cantor, Biochemistry, 1997, 36, 2339.3 D. P. Tieleman, S. J. Marrink and H. J. C. Berendsen, 8 R. D. Groot and P. B. Warren, J. Chem. Phys., 1997, Biochim. Biophys. Acta, 1997, 1331, 235. 107, 4423. 4 9 R. Goetz and R. Lipowsky, J. Chem. Phys., 1998, 108, 739 S. M. Willemsen, T. J. H. Vlugt, H. C. J. Hoefsloot and 7. 5 P. J. Hoogerbrugge and J. M. V. A. Koelman, 6 P. Español and P. Warren, Europhys. Lett., 1995, 30, 7 Y. Zhang, S. E. Feller, B. R. Brooks and R. W. Pastor, Europhys. Lett., 1992, 19, 155. 191. J. Chem. Phys., 1995, 103, 10252. B. Smit, J. Comp. Phys., 1998, 147, 507. 10 R. S. Cantor, Biophys. J., 1999, 76, 2625. Paper 9/06472I PhysChemComm © The Royal Society of Chemistry 1999

ISSN:1460-2733    

DOI:10.1039/a906472i

出版商:RSC    

年代:1999

数据来源: RSC