Langmuir 1996, 12, 5763-5767
5763
Computer Simulations of Spontaneous Vesicle Formation Ame´rico T. Bernardes† Institute for Theoretical Physics, Cologne University, D-50923, Ko¨ ln, Germany Received March 7, 1996. In Final Form: August 7, 1996X An improved version of an earlier Larson-type model of solution of amphiphiles is used to simulate the spontaneous formation of vesicles. Starting with a solution of amphiphiles either uniformly distributed or aggregated in an ordered way, self-organization into vesicles was observed by carrying out extensive Monte Carlo simulations. The aggregation of single-tailed amphiphiles and the role of the bending rigidity in the chains are also discussed.
1. Introduction Amphiphiles in aqueous solution aggregate in many ways, which will depend, for instance, on the temperature, pressure, concentration, and/or type of surfactant.1,2 Hydrophobic effect produces aggregates with the polar head outward in contact with the water molecules hiding the hydrocarbon chains. This hydrophobic effect has been explained as an entropic effect due to the structuring of the water around free surfactant molecules.3 Other interactions, like van der Waals, make no appreciable contribution to this effect.4 At low temperatures, amphiphiles tend to form continuous aggregates, among them membranes. The increasing interest in membranes comes from the fact that they form the basis of almost all life forms, as well as of technological applications.5 They can also be used as filters, simulating the action of the cell membranes. Biological membranes are generally bilayers formed by double chain amphiphiles (phospholipids) and stabilized by proteins. Artificial single or multilamellar membrane structures can be obtained with single or polychain amphiphiles.6 To avoid contact between the water and its edges, where the water molecules are in partial contact with the tails, the membranes may bend to form closed structures called vesicles or liposomes. Thus, vesicles separate different portions of the solution, and this property has important applications, such as drug delivery systems, the production of artificial cells (like artificial erythrocytes),7 or even for wastewater treatment.8 Vesicles are widespread and form spontaneously in nature; however, their artificial production requires the input of considerable mechanical energy and they are normally metastable.9,10 Recently, various methods were developed to obtain spontaneous formation of artificial † Present and permanent address: Departamento de Fı´sica, Universidade Federal de Ouro Preto, Campus do Morro do Cruzeiro, 34500-000 Ouro Preto/MG, Brazil: e-mail address,
[email protected]. X Abstract published in Advance ACS Abstracts, October 15, 1996.
(1) Tanford, C. The Hydrophobic Effect: Formation of Micelles and Biological Membranes; John Wiley: New York, 1980. (2) Gompper, G.; Schick, M. Self-Assembling amphiphilic systems; Academic Press: London, 1994. (3) Israelachvili, J. N.; Wennerstrom, N. J. Chem. Phys. 1992, 96, 520. (4) Cevc, G.; Marsh, D. Phospholipid Bilayers: Physical Principles and Models; John Wiley and Sons: New York, 1987. (5) Lipowsky, R. Nature 1991, 349, 475. (6) Kunitake, T.; Ando, R.; Ishikawa, Y. Mem. Kyushu Univ. 1986, 46, 221. (7) Gaber, B. P.; Schnur, J. M.; Chapman, D. Biotechnological Applications of Lipid Microstructures; Plenum Press: New York, 1988. (8) Sicchierolli, S. M.; Mamizuka, E. M.; Carmona-Ribeiro, A. M. Langmuir 1995, 11, 2991. (9) Kaler, E. W.; Kamalakara Murthy, A.; Rodriguez, B. E.; Zasadzinski, J. A. N. Science 1989, 245, 1371. (10) Cantu´, L.; Corti, M.; Musolino, M.; Salina, P. Europhys. Lett. 1990, 13, 561.
S0743-7463(96)00209-0 CCC: $12.00
vesicles. This was achieved by using mixtures of different types of amphiphiles9,11,12 or by using two-tailed amphiphiles with a large head.10 The same phenomenon was produced using “artificial” amphiphiles.13 Membranes are considered as two-dimensional flexible sheets and they can be modeled by the Helfrich bending Hamiltonian.14 In general the discussion of formation and change in shape of vesicles is based on the assumption that the membrane is already formed,15-20 thus not allowing the microscopic description of the aggregation mechanisms of amphiphiles forming vesicles. This important feature can be realized by the use of computer models. Computer simulations on lattice models have been used to study the aggregation of amphiphiles and related features such as the critical micelle concentration, transition temperatures, stability of aggregates, etc.21-29 (for a detailed review of Monte Carlo (MC) simulations on lattice models see ref 30). The micelles are the typical aggregate obtained in these computer simulations. Drouffe et al.31 observed the formation of vesicle structures in a molecular dynamics simulation, by using hard spheres representing the amphiphilic molecules. Brindle and Care32 have obtained a vesicle-like structure in a Larson-type model. In their work, heads and tails appear homogeneously (11) Kaler, E. W.; Herrington, K. L.; Kamalakara Murthy, A. J. Chem. Phys. 1992, 96, 6698. (12) Kondo, Y.; Uchiyama, H.; Yoshino, N.; Nishiyama, K.; Abe, M. Langmuir 1995, 11, 2380. (13) Kunitake, T. Angew. Chem. 1992, 104, 692. (14) Helfrich, W. Z. Naturforsch. C. 1973, 28, 693. For a detailed introductory description of this analytical treatment see: Peliti, L. Amphiphilic Membranes, lectures at 1994 Les Houches Summer School. (15) Safran, S. A.; Pincus, P. A.; Andelman, D.; MacKintosh, F. C. Phys. Rev. A 1991, 43, 1071. (16) Kumaran, V. J. Chem. Phys. 1993, 99, 5490. (17) Morse, D. C.; Milner, S. T. Europhys. Lett. 1994, 26, 565. (18) Seifert, U. Fluid membranes - Theory of vesicle conformations; Habilitation Theses, Institut Fu¨r Festko¨rperforschung, Ju¨lich, 1994. (19) Porte, G.; Ligoure, C. J. Chem. Phys. 1995, 102, 4290. (20) Fattal, D. R.; Andelman, D.; Ben-Shaul, A. Langmuir 1995, 11, 1154. (21) Larson, R. G.; Scriven, L. E.; Davis, H. T. J. Chem. Phys. 1985, 83, 2411. (22) Chowhury, D.; Stauffer, D. J. Chem. Phys. 1991, 95, 7764. (23) Larson, R. G. J. Chem. Phys. 1992, 96, 6904. (24) Jan, N.; Stauffer, D. J. Phys. I. 1994, 4, 345. (25) Bernardes, A. T.; Henriques, V. B.; Bisch, P. M. J. Chem. Phys. 1994, 101, 645. (26) Stauffer, D.; Woermann, D. J. Phys. II 1995, 5, 1. (27) Liverpool, T. B.; Bernardes, A. T. J. Phys. II 1995, 5, 1003. (28) Liverpool, T. B.; Bernardes, A. T. J. Phys. II 1995, 5, 1457. (29) Chowdhury, D.; Bernardes, A. T.; Stauffer, D. Int. J. Mod. Phys C, in press. (30) Liverpool, T. B. In Annual Reviews of Computational Physics; Stauffer, D., Ed.; World Scientific: Singapore, in press. (31) Drouffe, J. M.; Maggs, A. C.; Leibler, S. Science 1991, 254, 1353. (32) Brindle, D.; Care, C. M. J. Chem. Soc., Faraday Trans. 1992, 88, 2163.
© 1996 American Chemical Society
5764 Langmuir, Vol. 12, No. 24, 1996
Bernardes
distributed in the membrane, and it is not clear that this aggregate represents the equilibrium configuration. The case of stability of membranes was discussed by Liverpool and Bernardes27 using preformed membranes as the starting point of the simulations. More recently, the selforganization of vesicles was observed by using binary solutions of two-tailed amphiphiles initially uniformly distributed.33 The dynamical process of vesicle formation is found to be in agreement with the phenomenological description of aggregation.4 This process first starts with the formation of micelles, followed by aggregation of micelles which gives rise to prevesicles and membranelike structures and finally forming the vesicles. In this paper we discuss the spontaneous vesicle formation by using single or two-tailed molecules, with or without bending rigidity in the tails. We also performed simulations starting with a bilayer membrane. In all these cases vesicles turn out to be the most stable aggregate. 2. The Model In the present simulation the system is represented by a cubic lattice of L × L × Lz sites. An amphiphilic molecule is represented by la connected sites and a water molecule by a single site. Periodic boundary conditions are used in both x and y directions. To save computer time, the systems is confined by walls in the z direction. In this way, amphiphiles moving in the z direction cannot cross the walls, but this movement is allowed in the other two directions, namely, x and y (the directions x, y, and z here adopted are arbitrary, since the system is isotropic). To satisfy the excluded volume condition, each site may be occupied only once. In this paper, we have performed simulations with two types of amphiphiles: (1) Amphiphiles with one tail consist of a hydrophilic head (a “water-like” particle) with one site, a liaison site (that plays a noninteractive role), and a hydrophobic tail, with la - 2 sites. In this case we have used amphiphiles with la ) 6. (2) Amphiphiles with two tails consist of a head with two sites (defined in the middle of the polymer chain, namely, in the positions [la/2] and [la/2 + 1]), two liaison sites (beside the head towards the ends of the chain), and two tails with (la - 4)/2 sites each. For this case we use la ) 12. In this model a water particle is represented by a discrete variable Si ) +1 and the amphiphilic molecule by a string of connected sites (monomers) with values: +1, 0, -1, -1, -1, -1 (for the first case, one tail); and -1, -1, -1, -1, 0, +1, +1, 0, -1, -1, -1, -1 (for the second case, two tails). Our two-tailed amphiphile is similar to that “gemini” amphiphile described in the molecular dynamics simulations study of aggregation performed by Karaborni et al.34 In order to simulate the hydrophilic and hydrophobic effects, we assume that water-water and tailtail interactions are attractive and water-tail interaction is repulsive. So, the interaction energies are given by Eww ) Ett ) -Ewt ) - ( > 0). Because of the geometry of the lattice an amphiphile can only bend by 90°. In some simulations with two-tailed amphiphiles, we have included a bending stiffness in the tails and we have assumed that each bend costs an energy κB. Thus, the total energy of the system is given by Na
E)-
∑ SiSj + k)1 ∑ bkκB 〈ij〉
(33) Bernardes, A. T. J. Phys. II 1996, 6, 169. (34) Karaborni, S.; et al. Science 1994, 266, 254.
(1)
Figure 1. Two-dimensional representation of the amphiphiles movements used in the program: reptation, kink, and buckling. Connected circles mean amphiphile molecules with length la ) 6. Filled circles represent water. Note that all the movements are reversible.
where the sum 〈ij〉 in the interaction term is taken only over nearest neighbor sites. The second term represents the bending stiffness, where Na is the total number of amphiphiles in the solution and bk is the number of bends present in the tails of each amphiphile (in the present work we assume κB ) ). We define the temperature of the system in units of the energy interaction between sites t ) kBT/, where kB is the Boltzmann constant. Amphiphiles can move in the lattice changing position with water molecules through different types of movements: reptation, kink-jump, and buckling as is sketched in Figure 1 (in order to obey detailed balance a reverse buckling movement was included). A detailed description of these movements and technical discussion of the algorithm can be found in ref 28. In each attempt a neighboring site (or two sites, in the case of buckling movement) is chosen at random, and if it is occupied by a water molecule, one attempts to move the amphiphile to this new position in the following ways: (a) in the reptation, the chain slithers along its length; (b) in the kink-jump, a right angle on the amphiphile flips; (c) in the buckling two neighboring sites in the middle of the amphiphile spontaneously “wiggle” the amphiphile. In the reverse-buckling a wiggle is destroyed by pulling an amphiphile by its extremities. It should be noted that self-avoidance is automatically maintained by refusing all moves where the new position of any one of the sites on the chain was already occupied by an amphiphile site. Our Monte Carlo simulations are based on the Glauber algorithm. Different configurations of the system are generated according to a Boltzmann distribution probability at a given temperature.35 This is done by attempting to move an amphiphilic molecule. The change in energy of the system ∆E ) [Enew - Eold] is calculated and the move is accepted with probability p ) 1/(1 + exp(+∆E/kBT)). A Monte Carlo step here refers to one attempt to move all the amphiphiles, one at a time, according to their individual four types of possible movements. 3. Results In the present work, we performed extensive simulations to study the process of spontaneous vesicle formation. Firstly, we present the results obtained for a system with L ) 40 and Lz ) 50, with 533 two-tailed amphiphiles of la ) 12 (representing a fraction of φa ) 0.08 of the total lattice volume). The temperature was set to t ) 2.2, i.e., (35) Herrmann, D. W. Computer Simulation Methods in Theoretical Physics; Springer-Verlag: Berlin, 1990. Binder, K. In Monte Carlo Methods; Binder, K., Ed.; Springer-Verlag, Berlin, 1979.
Spontaneous Vesicle Formation
Langmuir, Vol. 12, No. 24, 1996 5765
Figure 2. 3D representation of the initial membrane (A) and of the vesicle (B, C, and D) obtained after 40 × 106 MC steps. The system size is 40 × 40 × 50, the temperature is t ) 2.2, and the amphiphile concentration is φ ) 0.08 (533 amphiphiles). Part B shows the entire vesicle; in part C the vesicle has been sectioned at the middle showing the right portion, and its inner side. The inner water is represented in part C by mid-gray circles. Finally, in part D only the middle layer is shown, where the membrane separating different domains can be better observed. Black circles represent the head particles, white ones represent the liaison particles, and light gray represent the tail particles. The z axis points upward in all figures.
below the condensation temperature tc (tc is the temperature below which the amphiphiles aggregate in just one cluster28). In this case, we did not take into account the bending energy; that means only the first term of eq 1 was used. In contrast with earlier simulations,33 we start the present one with the amphiphiles forming a membrane, as shown in Figure 2A. The vesicle obtained after 40 × 106 MC steps is shown in Figure 2B. In this figure, the water molecules surrounding the aggregate were omitted. Figure 2C shows a longitudinal section made in the middle of the vesicle. There, the water molecules into the vesicle
are shown as dark gray circles. The internal surface is basically formed by heads and liaison sites and encapsulates a volume of water which filled this internal portion. Figure 2D presents the middle layer in the longitudinal direction, where the vesicle membrane is clearly seen. As one can observe in these figures, the heads and liaison particles are mostly on the external and internal surfaces of the vesicle membrane, in contrast with Brindle and Care results.32 In their paper, the amphiphilic monomers are homogeneously distributed in the vesicle membrane. On the other hand, in the present work, the tails particles
5766 Langmuir, Vol. 12, No. 24, 1996
Figure 3. Chain concentration profiles obtained in layers parallel to the symmetry plane of the vesicle. The system is 40 layers long and has no walls in this direction. The concentrations of the different amphiphilic components as well as the total concentration are shown. The dip on the two upper curves represent the inner vesicle space filled by water molecules. The meaning of the symbols is provided in the legend.
are “enveloped” by the external surfaces. Due to the fact that this model is simulated in a cubic lattice with all the sites occupied by either water or amphiphilic particles and all the particles have the same “size”, it is unavoidable the contact between water molecules and tail monomers. In fact, the number of head and liaison particles cannot “cover” entirely the tails in this model. As we will discuss below, this configuration is more energetic than the membrane one. However, due to the entropy, this is the most favorable configuration for the temperature used in this simulation. The major length of this vesicle is oriented in the x direction, and the lower part is in contact with a wall (this is the reason for the “flatness” of the vesicle “lower” side). However in other simulations we have obtained different orientations as well as the vesicle in different regions of the system. This means that qualitatively the final results were not influenced by the restriction of the movement in the z direction. The thickness of the double-layer membrane (see Figure 2D) is close to the length of a tail (and not twice as usually described for phospholipids membranes). We have observed here the same kind of interpenetration of the tails as reported in the experiments of Cantu´ et al.10 when working with ganglioside GM3 molecules. If one looks at a cross section perpendicular to the major axis, one sees a symmetric profile, as shown in Figure 3. The concentrations of each amphiphile component were measured in layers parallel to the x-z plane (they are counted in the y direction). The total concentration of chain monomers is represented by a solid line. The dip on the middle represents the space filled with water molecules. This dip is also observed in the concentration of tail monomers, i.e., most of the tail monomers are in the “left-right” sides of the vesicle. The higher concentration of heads and liaison sites observed in the internal region is due to the fact that the internal surface is smaller than the external one and, as half of the heads and liaisons are in each surface, the concentration inside must to be higher than that outside. The explanation that the vesicles are formed through the membrane to minimize the interaction energy on the borders is not confirmed in our simulations. In fact, we always have a competition between the interaction energy and the entropy, and the less energetic configuration does not mean necessarily the best configuration. Figure 4 shows the total energy of the system versus time for three
Bernardes
Figure 4. Interaction energy for three different simulations. The solid and dotted lines represent two different simulations starting with amphiphiles uniformly distributed, while the dashed line represents the case where a membrane is the initial configuration.
cases: the first two configurations starting with amphiphiles uniformly distributed (with different seeds for the random number generator), and the third one starting with a membrane. The same system size at the same temperature, as well the same number of molecules, were used in these simulations. The three different systems converge to about the same value of final energy. As one can see, the membrane shape (represented by dashed line) is less energetic than the vesicle. Thus, for this temperature the vesicle is the most stable state. Presumably for lower temperatures, (when the entropic effect has less weight) membranes might be more stable, but we have not checked this assumption. Some possible ways for the dynamical process of vesicle formation have been described in the literature:4,5 membrane bending forming a closed surface, fusion of micelles with rearrangement of the amphiphilic molecules, and fusion of vesicles. As we have shown in a previous work, for an initial configuration of amphiphiles uniformly distributed in the solution, the process is basically characterized by the initial aggregation in small micelles. Then, when these micelles aggregate in larger clusters, a rearrangement is observed, leading to the formation of membrane-like structures (like those found by Cantu´ et al.10) or prevesicle aggregates (a small double-layer cluster with the inner heads compressed in the middle). Finally, when these aggregates come together, a new rearrangement takes place, leading to the formation of the vesicle. In the case starting with a membrane, the process is entirely different. Due to the fact that we are simulating a lattice model (and also a small system), it is not possible to observe the bending of the membrane surface, as usually described for large surfaces. What we have observed is a rearrangement of amphiphiles, with the size of the cluster (the number of molecules forming the cluster) being basically constant. The microscopic rearrangement of molecules found here is similar to that described by Drouffe et al.31 A simulation in a larger system (L ) Lz ) 50) with a higher amphiphilic concentration (φa ) 0.1, i.e., 1041 amphiphiles) for temperature t ) 2.2 showed that larger vesicles would be formed in the equilibrium. At 40 × 106 MC steps, two vesicles (with 630 and 321 amphiphiles) were formed but were still growing, as we observed in the dynamical evolution of the largest cluster. Thus, the largest vesicle reported above does not represent the largest that might be observed in an infinite system. The introduction of the bending energy slows down the system. However, it was again possible to obtain vesicles.
Spontaneous Vesicle Formation
Langmuir, Vol. 12, No. 24, 1996 5767
Figure 6. Dynamics of the aggregation process, showing the behavior of the largest cluster for different cases of two-tailed amphiphiles. Open and filled circles represent two different simulations starting with amphiphiles uniformly distributed; the solid line represents the case starting with a membrane and “plus” represent the case where the bending energy was assumed and amphiphiles were initially uniformly distributed. Figure 5. 3D representation of the vesicles obtained after 40 × 106 MC steps with 533 two-tailed amphiphiles, where bending energy was considered. The system size is 40 × 40 × 50, the temperature is t ) 2.2. Part A shows the two vesicles and part B shows only the middle layer. The symbols have the same meaning as those in Figure 2. In this case the vesicles were formed in the middle of the system, showing that the walls do not play any role in the vesicle formation. The z axis points upward.
Figure 5 shows the two small vesicles obtained after 40 × 106 MC steps for the 40 × 40 × 50 system with amphiphiles la ) 12 long, by taking into account the bending energy. The simulation started with amphiphiles uniformly distributed. Figure 5A shows the two globular vesicles and Figure 5B shows the middle layer of these aggregates. The bending energy forces the external surface to be flatter than that without bending energy. Small vesicles obtained where the bending energy was not assumed were more “spherical”. This happens because tails tend to “prefer” a straight conformation. The bending energy also increases the rigidity of the aggregate which leads to the slowdown of the system. In fact, as we have observed in a previous work, the clusters perform an “amoeba-like” movement.33 In order to perform this movement, sometimes an amphiphile must to “try” to leave the cluster by “pulling” it. As the bending energy forces the tails to be more straight, this makes different types of amphiphilic movements and therefore the movement of the whole cluster more difficult. Figure 6 shows the different pathways describing the dynamical evolution of the largest cluster for the different simulations of two-tailed amphiphiles described above. For fixed time intervals (here 400 × 103 MC steps) we count all the amphiphile clusters. An amphiphile belongs to a cluster if any of its monomers is a neighboring site to any of the monomers of amphiphiles already counted as part of the cluster. Circles (open and filled) represent amphiphiles without bending energy and initially uniformly distributed. The membrane initial configuration case is represented by a solid line (notice that the largest cluster contains all the amphiphiles and keeps all the time about the same size). The “plus” represents the case where bending energy was assumed. As one can see, the introduction of bending energy slows down the vesicle formation. Finally, we also performed simulations for single-tailed amphiphiles. For the same system size we have used the same amphiphilic concentration φ ) 0.08 (1066 amphiphiles) at the same temperature t ) 2.2. After 40 ×
106 MC steps two vesicles were formed (with 286 and 776 molecules each). The two vesicles interchange molecules and presumably they will form only one vesicle at a later time. The aggregation process with two-tailed amphiphiles is faster because they can “pull” faster the whole cluster. The buckling move may play an important role in this process, due to the fact that in the single-tailed case it simply bends the amphiphile inner portion. In contrast with this result, when studying the stability of bilayers in a binary mixture, Liverpool and Bernardes27 observed that below tc ) 2.3 the bilayer was stable, but above tc the bilayer dissolved. They used single-tailed amphiphiles in their simulations with la ) 6 and monomers strength values +2, +1, 0, -1, -1, -1 on systems of 48 × 48 × 52 for times up to 106 MC steps. In that study the simulations always started with a membrane. However, when the amphiphiles were uniformly distributed at the beginning of the simulations, they aggregated in a clumsy cluster. As one can observe, they used amphiphiles with a strong head-head interaction, in contrast with those used in this paper or in the Brindle and Care work.32 Moreover, they did the simulations by using just the reptation movement. Presumably, the strong interaction between heads and the using of just reptation did not allow the system to “explore” different ways of aggregation. Simulations with shorter two-tailed amphiphiles la ) 10 has also been performed, yielding vesicle formation. It is interesting to note that we have not observed different dynamical behaviors in the aggregation process, such as reported by the Shell group,34 when comparing the results obtained with MD simulations for single-tailed and twotailed amphiphiles. Basically, we observed the slowing down in the vesicle formation with single-tailed, but no such “tree-like” micelles have been formed. Each of the reported simulations spent around 5-6 days of CPU in a 2GB DEC Alpha AXP/7000 of the RRZ/Cologne University. The complete Fortran code is available upon request. Acknowledgment. I would like to thank D. Chowdhury, D. Stauffer, T. B. Liverpool, and C. Shida for many stimulating discussions and suggestions and O. F. de Alcantara Bonfim for carefully reading the manuscript. This work is partially supported by the Brazilian Agency Conselho Nacional para o Desenvolvimento Cientı´fico e Tecnolo´gicosCNPq. LA960209J
