Spontaneous Vesicle Self-Assembly: A Mesoscopic View of

ARTICLE pubs.acs.org/Langmuir

Spontaneous Vesicle Self-Assembly: A Mesoscopic View of Membrane Dynamics Julian C Shillcock* MEMPHYS, University of Southern Denmark, Campusvej 55, 5230 Odense M, Denmark ABSTRACT: Amphiphilic vesicles are ubiquitous in living cells and industrially interesting as drug delivery vehicles. Vesicle self-assembly proceeds rapidly from nanometer to micrometer length scales and is too fast to image experimentally but too slow for molecular dynamics simulations. Here, we use parallel dissipative particle dynamics (DPD) to follow spontaneous vesicle self-assembly for up to 445 μs with near-molecular resolution. The mean mass and radius of gyration of growing amphiphilic clusters obey power laws with exponents of 0.85 ( 0.03 and 0.41 ( 0.02, respectively. We show that DPD provides a computational window onto fluid dynamics on scales unreachable by other explicit-solvent simulations.

’ INTRODUCTION Phospholipid bilayer membranes are fundamental for the stability and function of living cells.1 They protect the cell from the external world while still allowing it to communicate, expel waste, and take in nutrients. Artificial polymeric membranes are becoming increasingly important for industrial applications ranging from drug delivery vehicles2 to fuel cell membranes.3 Phospholipids in bulk solution spontaneously form complex supramolecular aggregates including micelles, bilayers, and vesicles driven by the amphiphilic nature of the constituent molecules, in which a hydrophilic segment is covalently bonded to a hydrophobic segment. Because the fast kinetics of the early stages of such supramolecular rearrangements are experimentally invisible, computer simulations such as molecular dynamics4 (MD) have been used to visualize them. However, the complexity of the interatomic force fields limits them to membrane patches5 and vesicles6 some tens of nanometers in extent for times of tens of microseconds or less. Amphiphile self-assembly on length scales beyond 20 nm is a computationally slow process that is currently only possible using coarse-grained models. Various coarse-grained simulation techniques have been developed for phospholipid membranes7 that allow them to reach length and time scales inaccessible to atomistic MD. They achieve this feat largely by grouping atoms into particles, thereby reducing the number of degrees of freedom that must be integrated in the equations of motion. These coarse-grained particles interact by softer effective forces than the typical interatomic force fields used in MD, which provides another speed up. The benefits of large-scale simulations of amphiphilic membranes are, however, more significant than simply reducing the simulation time. Fluctuation-induced interactions between distant parts of a membrane coupled with hydrodynamic modes in the r 2011 American Chemical Society

solvent give rise to dynamical correlations that can guide a large system along a different kinetic pathway from that followed by a small system. Dissipative particle dynamics (DPD) is a coarsegrained technique invented almost 20 years ago8,9 to study hydrodynamic phenomena in fluids that require explicit solvent at a reduced computational cost than using MD. It is characterized by the use of soft, short-ranged interparticle forces and a thermostat that conserves momentum locally. It has been used to study the microphase separation of polymer melts,11 phospholipid/water phase diagrams,12 the self-assembly and properties of lipid bilayers,13,14 diblock copolymer membranes,15 vesicle fusion,16
19 budding of two -component vesicles,20 and many more. However, published results of DPD simulations of amphiphilic membranes (see refs 21 and 22 for two recent extensive reviews) are limited to systems less than 50 nm in spatial extent except for one or two exceptions.23,24 In order to go beyond this size, within a reasonable computational time frame, a parallel implementation of DPD is required. Note that the term DPD is sometimes used to refer only to the Galilean-invariant thermostat, but we use it here to mean the thermostat and the set of soft forces defined in the original scheme.8
10 Atomistic MD has also been used to study the formation of small surfactant micelles,25 and 10
20 nm diameter vesicles,6,26 and coarse-grained MD has been used to follow the fusion of such vesicles.27,28 However, formation of multiple vesicles on the scales we observe has only been studied before using dynamic self-consistent field theory,29 which does not include hydrodynamic effects that may influence the kinetics of structure evolution in the system. Also, the high curvature of small vesicles Received: April 28, 2011 Revised: November 2, 2011 Published: November 22, 2011 541

dx.doi.org/10.1021/la2033803 | Langmuir 2012, 28, 541–547

Langmuir

ARTICLE

Figure 1. Time series of snapshots of vesicle self-assembly from an initially random dispersion of 58 240 amphiphiles in 5 000 000 water particles contained in a (90 nm)3 simulation box (80 mM amphiphile concentration). Amphiphiles have the linear architecture H3T6, and the hydrophilic H particles are colored red and the hydrophobic tail particles orange except that the terminal tail particle is green. These snapshots illustrate the three stages of growth. Initially, the amphiphiles form small micelles that diffuse, merge, and transform into planar bilayer patches (a, log(time) = 9.9, ∼2 μs). Next, a middle phase of growth appears during which the bilayer patches grow and curl up into vesicles that subsequently continue to grow by fusion with each other and the remaining bilayer sheets (b, log(time) = 13.6, ∼80 μs). Eventually the system is dominated by closed vesicles that diffuse around but only occasionally fuse as the vesicle fusion process requires times longer than the current simulations (c, log(time) = 14.3, ∼162 μs). The duration of the stages depends on the amphiphile concentration as discussed in the text.

makes them unsuitable models for vesicles in biology (typically 40 nm or larger in diameter) and the pharmaceutical industry (drug delivery vesicles are hundreds of nanometers in diameter). We observe that vesicle self-assembly passes through three successive stages, which are dominated by aggregates that grow by distinct kinetic mechanisms. Initially, and very rapidly, small micelles form and merge into larger micelles. Next, the micelles transform into or merge with fluctuating, quasi-planar bilayer sheets that curl up into closed vesicles. During this phase we extract the growth laws for the mean cluster mass and radius of gyration. Once the system is dominated by closed vesicles, the size distribution reaches a plateau that persists for the remainder of our simulations, as the fusion of large vesicles takes place on a still longer time scale. The second phase of growth may be relevant to the loading of vesicular drug delivery vehicles, which depends on getting the largest amount of active material inside the vesicles as they form.30 Understanding the kinetics of vesicle formation during this second phase of growth is necessary before one can rationally optimize this process.

sodium dodecylphosphate. Higher concentrations lead to faster vesicle formation and so reduce the simulation time required. In DPD, as in coarse-grained MD, the number of degrees of freedom in each molecule is reduced by combining several atoms or atomic groups into particles that interact via an effective force field. Several water molecules are grouped into one water particle, and the interactions are chosen so as to reproduce the compressibility of water at room temperature.10 Several methyl groups in a hydrocarbon chain are grouped into a single hydrophobic particle. We do not try to capture the chemical details of a particular lipid or polymer, for which DPD is not a suitable technique, but choose a molecular architecture that has been shown in previous work to lead to stable bilayer structures.14 It is representative of a typical amphiphile and has a linear architecture H3T6, in which three hydrophilic “Head” particless (H) are attached to a chain of six hydrophobic “Tail” particles (T). All particles have the same radius a0, which is the range of the nonbonded forces. To reduce the time required for force calculations, no bending stiffness potential14 is imposed on the amphiphiles beyond that which originates in the nonbonded interactions of the particles. Water is represented as a single particle (W). Once the molecular architecture, concentration, and force field have been specified, the simulation evolves by integrating Newton’s laws of motion for all particles. All simulations are carried out at constant temperature (300 K). Because DPD is a coarse-grained technique, it is necessary to assign mass, length, and time scales to the simulated quantities by comparing suitable observables to experimental values. As in previous work,14,19

’ RESULTS AND DISCUSSION We performed parallel DPD simulations of vesicle selfassembly from an initially random dispersion of 30 208 and 58 240 model amphiphiles and more than 5 000 000 water particles in a (90 nm)3 box. These represent concentrations of 41 and 80 mM, respectively, which are well above the micro- to millimolar CMC of single-tailed amphiphiles such as lysolipids or 542

dx.doi.org/10.1021/la2033803 |Langmuir 2012, 28, 541–547

Langmuir

ARTICLE

we assume that all particles have equal mass and fix the length and time scales by simulating a tensionless planar membrane composed of H3T6 amphiphiles and comparing their area per molecule and inplane diffusion constant with those of the typical phospholipid dimyristoylphosphatidylcholine. This yields the values a0 = 0.7 nm and 1 time step is 0.1 ns. The vesicle formation simulations then represent between 172 and 445 μs of real time. More details about DPD simulations of amphiphiles are given in the literature.10,14,21 Vesicle self-assembly is observed to pass through three distinct stages that are independent of the amphiphile concentrations we studied, although the duration of each stage varies with concentration. Starting from a random dispersion of the amphiphiles in bulk water, there is a period, which lasts up to a few microseconds (Figure 1a), that results in formation of many small, spherical micelles that gradually merge and transform into quasi-planar bilayer patches. We note here that in earlier atomistic MD simulations31 54 surfactants initially dispersed in a simulation box with linear dimension 9 nm formed a single micelle within a few nanoseconds. Small micelles also form as rapidly in our simulations, but it is the duration over which they remain micelles that characterizes this stage of the vesicle formation process. This duration is controlled by the time required for micelles to grow large enough to transform into planar membrane patches. Second, there is a longer period during which the micelles and planar bilayer patches diffuse around and slowly merge into larger bilayer patches that curl up into closed vesicles (Figure 1b). This stage lasts from 50 to 400 μs depending on the concentration and is shorter for higher concentrations. It encapsulates the potentially important process of vesicle formation and closure. We note that in the concentration range we studied vesicles and planar bilayer patches coexist for almost all of this stage. The vesicles appear to be unilamellar and small at lower concentrations but larger at higher concentrations. Finally, the system is dominated by closed vesicles (Figure 1c), whose diffusion and merging take place on time scales longer than our simulations. Typically, we observe that about 80
90 vesicles and membrane fragments remain in a box of (90 nm)3 and 20
30 in a box of (70 nm)3, and these numbers are broadly independent of the amphiphile concentration. These are unlikely to be the final equilibrium states, however, as fusion of vesicles at these concentrations appears to take place on longer times than we simulated. This regime could also be probed using parallel DPD if a larger number of processors were dedicated to the task. In order to quantify the vesicle formation process, we measure geometrical properties of the growing clusters. As the initially dispersed amphiphiles aggregate they form clusters (micelles, bilayer patches, etc.) that are surrounded by an expanding region of water. This is seen in Figure 1a
c, where the solvent-filled spaces between clusters grow larger as time passes. We assign all the amphiphiles to a set of discrete clusters labeled with the integers, 1, 2, 3, etc., as follows. Two amphiphiles belong to the same cluster if their first head particles are separated by less than a fixed distance Dmax/a0 = 2. This value of Dmax was chosen because smaller values around Dmax/a0 ≈ 1 cause the algorithm to treat thermally fluctuating amphiphiles as (transiently) leaving the cluster and leads to single clusters being counted as distinct clusters. Likewise, larger values of Dmax/a0 ≈ 5 result in discrete clusters separated by a thin water gap being counted as one. We verified that the mean quantitative properties of the clusters are unchanged if we vary Dmax/a0 around 2 and if we use the first head particle or the last tail particle to define the location of the amphiphiles. We iterate over all the amphiphiles, calculating the

Figure 2. (a) Comparison of the mean mass of amphiphilic clusters as a function of time with that of oil droplets aggregating in water. Two independent simulations with amphiphile concentrations of 41 (lower two curves) and 80 mM (middle two curves) are shown, which correspond to 30 208 and 58 240 amphiphiles. There are 37 496 oil molecules (upper curve) corresponding to 10 vol %. The best-fit straight lines to the oil droplet and vesicle growth curves are shown. An explanation of which points are used to calculate the growth exponents and the errors is given in the text. The best-fit exponents for the curves shown are 1.0 for the oil and 0.88 for the 80 mM amphiphile concentrations and 0.87 for the 41 mM amphiphile concentrations, as shown by the straight lines. Note that these linear regression fits are only to the curves shown. The results from the complete set of 4 simulations at each amphiphile concentration are described in the text. The vesicle growth law is independent of the amphiphile concentration in the range studied. The amphiphilic clusters grow more slowly than oil droplets with a growth exponent that is close to the value of 6/7 ∼ 0.857 predicted by the heuristic model in the text. Note that the final state of the vesicles is not an equilibrium state as vesicle fusion takes place on longer times than we simulated. (b) Comparison of the mean radius of gyration of the amphiphilic clusters and oil droplets for the same systems as in Figure 2a. The growth exponent for the radius of gyration is 0.33, as expected qfor the oil droplets, while it is 0.42 for both the 80 and 41 mM amphiphile concentration systems, which is close to the predicted value of 3/7 ∼ 0.429.

distances from each one to all preceding ones, and if this distance lies within Dmax of an amphiphile within an existing cluster or another single amphiphile we add it to the cluster or construct a new cluster, respectively. As clusters are constructed, it sometimes happens that two clusters are created that are found to be connected by an amphiphile checked later in the iteration loop, in 543

dx.doi.org/10.1021/la2033803 |Langmuir 2012, 28, 541–547

Langmuir

ARTICLE

which case we transfer all the amphiphiles in the higher numbered cluster to the lower numbered one. Once all molecules have been assigned to clusters, we relabel the clusters with contiguous integers and analyze their properties. We count the total number of clusters (micelles, bilayer patches, and vesicles) at fixed time intervals and measure their mean mass (i.e., number of amphiphiles per cluster) and radius of gyration. We use periodic boundary conditions in the simulations, which means that particles that leave the simulation box through a face enter it immediately at the opposite face. Fragments of a cluster that span the periodic boundaries are counted as one cluster in the analysis. We predict the time dependence of cluster growth by a heuristic argument that accounts for diffusion and coalescence of growing clusters. Let there be N amphiphilic clusters in a volume of linear extent L (in units of the particle radius a0), the mean volume per cluster is then L3/N, and the mean distance between the centers of mass of any two nearest-neighbor clusters is L/N1/3. Assuming that a cluster of radius R has a diffusion constant D given by the Stokes
Einstein relation D = kBT/(6πηR), where T is the temperature, kB is Boltzmann’s constant, and η is the viscosity of water, the mean time, t, for one cluster to diffuse to and merge with another cluster is given by ÆL2/N2/3æ ≈ Dt. Substituting for D gives t ≈ R/N2/3. Note that we assume that the time required for two clusters to merge once contact has been made is negligible compared to the time required for them to diffuse to each other. This assumption is seen to be true by inspection of successive snapshots in the simulations. The total number of amphiphiles in all clusters is constant, and if we assume that the system contains only quasi-planar bilayer patches and closed vesicles, for which the area grows as the square of their linear size, we have the relation N ≈ R
2. Therefore, the mean cluster radius grows with time as R(t) ≈ t3/7 and their mean mass as M(t) ≈ t6/7. We measure the cluster radius of gyration, Rg, in the simulations, as it is a more general quantity for arbitrary shapes than radius, which tends to presuppose a spherical shape. Their scaling with time is expected to be the same. A similar argument can be applied to the coalescence of oil droplets in bulk water except that the constraint of a constant number of oil molecules distributed among N droplets leads to N ≈ R
3, leading to R(t) ≈ t1/3 and M(t) ≈ t. In order to check this, we simulated a mixture of oil molecules composed of a single particle (T) with the same interactions as the amphiphile tail particles and water (W) at a 10% number fraction of oil and measured the exponents for the mean oil droplet mass and radius of gyration. The data shown in Figure 2a and 2b for oil droplets in water are from a simulation of 37 496 oil particles and 337 496 water particles in a simulation box (35 nm)3. The mean oil droplet mass and radius of gyration closely follow the predicted curves except at late times. The flattening of the growth curves seen in this limit appears because there are only two large oil droplets remaining that diffuse too slowly to merge during the simulation and whose masses are therefore constant for a long time. The tiny fluctuations to smaller mass at late times result from oil particles detaching from a droplet and diffusing around in the solvent before rejoining a droplet. The exponents obtained for oil droplet growth are the same as those predicted by Lifshitz
Slyozov theory32 for the ripening of droplets of one phase growing in a second phase via evaporation and condensation of single molecules. This process may be active during oil droplet growth in our simulations but cannot be dominant in the vesicle selfassembly process as it would lead to larger growth exponents for vesicle mass and radius of gyration than we observe.

Figure 3. Growth of the mean mass of amphiphilic clusters as a function of time for amphiphilic systems in boxes with volume (114 nm)3 and (140 nm)3. The numbers of amphiphiles in each system are (starting with the lowest curve) as follows: 58 368 in (140 nm)3 (21 mM), 170 000 in (140 nm)3 (62 mM), 86 528 in (112 nm)3 (62 mM), 180 224 in (140 nm)3 (66 mM), and 102 912 in (112 nm)3 (73 mM). As a guide to the eye, two straight lines with slope 0.82 are overlaid on the lowest and highest curves (21 and 73 mM, respectively) to show that all systems obey the same growth law to the accuracy of the simulations. The vesicle simulations, especially the most dilute system (21 mM), deviate from the expected growth law at early times because the rapidly formed micelles take longer to form the bilayer patches and partial vesicles that are necessary for the middle stage of growth as described in the text. Note that the deviation between the growth curves at a fixed concentration is less than the symbol size.

Figure 2a and 2b shows the time dependence of the mean mass and mean radius of gyration of the growing amphiphilic clusters compared to that of oil droplets in bulk water. The slopes of the linear parts of the curves in Figure 2a and 2b give the growth exponents for the mean cluster mass, α, and radius of gyration, β, respectively, from MðtÞ ≈ t α R g ðtÞ≈t β with α = 0.85 ( 0.03 and β = 0.41 ( 0.02 for amphiphiles and α = 1.02 and β = 0.34 for oil droplets, respectively. The errors quoted are estimated from 4 independent simulations at each of the two amphiphile concentrations. We do not calculate standard deviations for the oil droplet exponents as we only have 2 data points and these results are intended for comparison only. Also note that only 2 vesicle systems at each concentration and 1 oil droplet are shown in the figure for clarity. All of the exponents extracted from the complete set of 13 vesicle simulations (four of which are shown in Figure 2a and 2b, four of which are included in the exponent averages but not in the figure, and five in Figure 3) and two oil droplet simulations, together with the system size and run length of the simulations are summarized in Table 1. The growth exponents for the oil droplets are calculated from the first 500 000 time steps (log(time) < 13.1) to exclude the plateau at large times when only one or two large droplets are present. For the vesicles at 80 mM concentration, we use the range 50 000
1 000 000 time steps (10.8 < log(time) < 13.8), while for vesicles at 41 mM concentration we use the range 100 000
2 000 000 (11.5 < log(time) < 14.5) to exclude the 544

dx.doi.org/10.1021/la2033803 |Langmuir 2012, 28, 541–547

Langmuir

ARTICLE

Table 1. Summary of the Simulations Performed To Extract the Growth Exponents for Amphiphilic Vesicles and Oil Droplets in Watera amphiphile number

box size (nm)

concentration (mM)

α

β

run time (106 steps)

30 208

90

41

0.885

0.433

2.2

30 208 14 400

90 70

41 42

0.840 0.789

0.412 0.390

2.18 3.2

14 400

70

42

0.849

0.424

4.45

58 240

90

80

0.878

0.432

1.72

58 240

90

80

0.876

0.420

1.73

27 776

70

81

0.855

0.374

3.18

27 776

70

81

0.845

0.400

4.33

58 368

140

21

0.812

0.406

3.06

170 000 86 528

140 112

62 62

0.877 0.871

0.424 0.432

3.06 3.7

180 224

140

66

0.827

0.408

1.6

102 912

112

73

0.846

0.421

1.45

37 496

35

874

0.998

0.326

1.0

37 496

35

874

1.036

0.356

2.0

The number of amphiphiles (or oil molecules), box size, and concentration are listed in the first three columns. The growth exponents for the mean cluster mass (α) and mean cluster radius of gyration (β) are given in the fourth and fifth column, and the last column shows the total number of time steps in each of the runs. The first eight rows show the vesicle data, only four curves of which are shown in Figure 2a and 2b, while the next five rows are data for the largest vesicle systems, all of which are shown in Figure 3. The corresponding results for the oil droplets growing in water are shown in the final two rows of the table. All simulations were performed at constant particle number, volume, and temperature. a

times, we do not see the larger systems enter the final stage of growth by vesicle fusion as this requires longer times that we have been able to simulate. The runs in Figure 3 require between 4 and 7 days of computing time on 512 processors, corresponding to 6
10 CPU-years each. The exponents for the oil droplets agree very well with those expected from the coalescence of droplets and classical Lifshitz
Slyozov theory,32 which are α = 1 and β = 1/3. The growth exponents for the mass and radius of gyration of the amphiphilic clusters also agree well with those expected from the above heuristic argument (6/7 ∼ 0.857, 3/7 ∼ 0.429). This indicates that the vesicles grow by coalescence of bilayer fragments (not necessarily planar) and not by the evaporation/condensation mechanism of Lifshitz
Slyozov theory. Visual inspection of the simulation snapshots shows that the middle phase of growth consists of vesicles or partially formed vesicles together with fluctuating, quasi-planar bilayer patches (see Figure 1b). The duration of this middle phase of vesicle formation depends on the total amphiphile concentration. It ranges from 50 μs to more than 400 μs, being shorter for higher concentrations as the vesicles more rapidly reach the stage in which they are large and widely separated in space. Subsequent growth of the vesicles can only occur by their diffusing to each other and fusing, which appears to take place on longer times than we probed here. However, our results show that the growth exponent for vesicle formation is independent of the amphiphile concentration in the range studied. We also checked that the mean mass of the growing clusters in the middle stage is linearly proportional to the square of their radius of gyration (data not shown), showing that the vesicles are approximately spherical and not cylindrical. It is perhaps worth making a comment about why parallel DPD is particularly suitable for simulating complex fluids or soft matter systems. As computing technology tends toward increasing the number of cores per processor chip instead of simply increasing the processor clock frequency, it is important that

initial micelle-dominated state and the plateau at large times when the system is dominated by closed vesicles that fuse on time scales longer than our simulations. As a check that our results are not limited only to small system sizes, we also simulated vesicle formation in larger systems. Figure 3 shows the increase in the mean vesicle mass as a function of time for systems containing from 58 368 amphiphiles in (140 nm)3 (which corresponds to 21 mM) to 180 224 amphiphiles in (140 nm)3 (which corresponds to 66 mM). We overlaid two straight lines with slope 0.82 on the upper and lower curves in Figure 3 to indicate that the growth law for the larger systems is independent of concentration and simulation box size and, although slightly smaller, within one standard deviation of that obtained for the (90 nm)3 systems. Inspection of Figures 2a, 2b, and 3 shows that the mean cluster mass and radius of gyration of the amphiphilic clusters grow more slowly at early times than the above argument would suggest and that this deviation is not seen for the oil droplets. We hypothesize that it is due to the distinct mechanisms by which micelles grow compared to oil droplets. The free energy of a system of oil droplets decreases continuously as the oil phase separates from the water. However, the free energy of a micelle forming out of a random dispersion of amphiphiles reaches a local minimum at some aggregation number as it becomes harder to pack extra amphiphiles into an existing quasi-spherical micelle. Once the system has established this micelle size distribution, further growth occurs mainly by micelles diffusing across the solvent gaps separating them and fusing with each other, a process that transiently perturbs the micelle’s preferred packing. The micelle-dominated stage is therefore metastable, and the mean cluster mass grows more slowly than in the secondary stage of vesicles and bilayer patches. This hypothesis is supported by the data in Figure 2b, which shows the radius of gyration of the amphiphilic clusters is quite shallow at early times, which is consistent with a metastable state of spherical micelles. At late 545

dx.doi.org/10.1021/la2033803 |Langmuir 2012, 28, 541–547

Langmuir computer simulation codes take maximum advantage of this trend. Use of Infiniband-connected commodity processors combined with the Message Passing Interface protocol to handle the interprocessor communication provides terascale computing power at a fraction of the cost of a dedicated supercomputer. In the absence of electrostatic forces, the nonbonded forces in DPD are short ranged. A large simulation box can be split up into a cuboidal grid of small volumes, each of which is assigned to a separate processor, a procedure known as spatial domain decomposition. All particles in the volume of space allocated to a processor are “owned” by that processor. The short range of the soft DPD forces ensures that only particles located within one force cutoff distance of a processor’s boundary can interact with particles owned by neighboring processors. All other particles within a processor’s space only interact with particles owned by the same processor. Crucially, if the linear size of each processor’s space is at least 20 times the force cutoff distance, the cost of the interprocessor communication is almost negligible compared to the force calculations each processor has to perform for particles in its own volume. We have shown that parallel DPD on a modest cluster of 64 processors can simulate the spontaneous selfassembly of vesicles from 58 240 amphiphiles in a box of size (90 nm)3 in less than 1 week of real time. Larger systems of 180 224 amphiphiles in a (140 nm)3 exhibit the same behavior but require 1 week on 512 processors, corresponding to 10 CPUyears per simulation. We observe a linear speedup of the simulation code up to 512 processors and expect this speed up to continue indefinitely. This advantage is not unique to DPD and also applies to parallel molecular dynamics, for example. However, the force fields used in MD are computationally expensive, which means that the volume assigned to each processor must be smaller to allow the simulation to proceed in a reasonable (real) time, and this reduces the ratio of the local force calculation to messaging time resulting in a lower gain for each added processor.

ARTICLE

Yet the simple question of whether nanometer-sized, hydrophobic aggregates do or do not self-assemble in the hydrophobic interior of a lipid membrane has not yet been definitively answered. Previous MD simulations36 indicate that fullerenes do not aggregate inside model membranes, but the behavior of larger aggregates is unknown. The relevant length and time scales can easily be reached using parallel DPD.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: Julian.Shillcock@epfl.ch. Present Addresses

Blue Brain Project, EPFL, Quarti
er de l’innovation, CH-1015 Lausanne, Switzerland.

’ ACKNOWLEDGMENT The author would like to thank John Ipsen, Chris Lagerholm, Mohamed Laradji, and Michael Lomholt for very helpful discussions during this work and gratefully acknowledges support from MEMPHYS and use of the Danish Centre for Scientific Computing’s “Horsehoe” supercluster at the University of Southern Denmark. MEMPHYS is supported by the Danish National Science Foundation. ’ REFERENCES (1) Alberts, B.; Bray, D.; Lewis, J.; Raff, M.; Roberts, K.; Watson, J. D. Molecular Biology of the Cell, 2nd ed.; Garland Publishing: New York, 1989. (2) Meng, F.; Engbers, G. H. M.; Feijen, J. J. Controlled Release 2005, 101, 187–198. (3) Dorenbos, G.; Suga, Y. J. Membr. Sci. 2009, 330, 5–20. (4) Rapaport, D. C. The Art of Molecular Dynamics Simulation; Cambridge University Press: Cambridge, 1995. (5) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. J. Phys. Chem. B 2007, 111, 7812–7824. (6) Marrink, S. J.; Mark, A. E. J. Am. Chem. Soc. 2003, 125, 15233–15242. (7) Bennun, S. V.; Hoopes, M. I.; Xing, C.; Faller, R. Chem. Phys. Lipids 2009, 159, 59–66. (8) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Europhys. Lett. 1992, 19, 155–160. (9) Espag~ nol, P.; Warren, P. B. Europhys. Lett. 1995, 30, 191–196. (10) Groot, R. D.; Warren, P. B. J. Chem. Phys. 1997, 107, 4423–4435. (11) Groot, R. D.; Madden, T. J.; Tildesley, D. J. J. Chem. Phys. 1999, 110, 9739–9749. (12) Kranenburg, M.; Smit, B. J. Phys. Chem. B 2005, 109, 6553–6563. (13) Venturoli, M.; Smit, B. Phys. Chem. Commun. 1999, 2, 45–49. (14) Shillcock, J. C.; Lipowsky, R. J. Chem. Phys. 2002, 117, 5048–5061. (15) Ortiz, V.; Nielsen, S. O.; Discher, D. E.; Klein, M. L.; Lipowsky, R.; Shillcock, J. C. J. Phys. Chem. B 2005, 109, 17708–17714. (16) Shillcock, J. C.; Lipowsky, R. Nat. Mater. 2005, 4, 225–228. (17) Liu, Y.-T.; Zhao, Y.; Liu, H.; Liu, Y.-H.; Lu, Z.-Y. J. Phys. Chem. B 2009, 113, 15256–15262. (18) Wu, S.; Guo, H. J. Phys. Chem. B 2009, 113, 589–591. (19) Grafm€uller, A.; Shillcock, J. C.; Lipowsky, R. Biophys. J. 2009, 96, 2658–2675. (20) Laradji, M.; Kumar, P. B. S. J. Chem. Phys. 2005, 123, 224902. (21) Venturoli, M.; Sperotto, M. M.; Kranenburg, M.; Smit, B. Phys. Rep. 2006, 437, 1–54.

’ CONCLUSIONS In summary, parallel dissipative particle dynamics simulations of the spontaneous self-assembly of amphiphiles into vesicles reveals a stage in which vesicles coexist with fluctuating, planar bilayer patches. The mean mass and radius of gyration of these vesicles/patches grow with power laws distinct from those found for oil droplets coalescing in water. The measured exponents of 0.85 ( 0.03 and 0.41 ( 0.02 agree well with a simple model that predicts 6/7 and 3/7, respectively. This stage may be important for optimizing the encapsulation efficiency of vesicular drug delivery vehicles30 but is currently inaccessible to experimental techniques. For example, fluorescence correlation spectroscopy33 with a sampling frequency of 100 kHZ does not yield enough data points during the few hundred microsecond duration of the stage to construct the correlation functions required to measure the vesicle size distribution. The simulations therefore predict the vesicle growth dynamics in a regime complementary to experiments. Parallel DPD can also be applied to other mesoscopic membrane processes. Its ability to probe length and time scales orders of magnitude larger than other explicit-solvent, particle-based techniques allows entropic forces that are negligible in smaller systems to be observed. The increasing commercial use of carbon nanotubes and spherical fullerenes (C60) in many materials has made the problem of understanding their interactions with biological tissue, in particular, cellular membranes, acute.34,35 546

dx.doi.org/10.1021/la2033803 |Langmuir 2012, 28, 541–547

Langmuir

ARTICLE

(22) Marrink, S. J.; de Vries, A. H.; Tieleman, D. P. Biochim. Biophys. Acta 2009, 1788, 149–168. (23) Shillcock, J. C.; Lipowsky, R. J. Phys.: Condens. Matter 2006, 18, S1191–S1219. (24) F€uchslin, R. M.; Maeke, T.; McCaskill, J. S. Eur. Phys. J. E 2009, 29, 431–448. (25) Jorge, M. Langmuir 2008, 24, 5714–5725. (26) de Vries, A. H.; Mark, A. E.; Marrink, S. J. J. Am. Chem. Soc. 2004, 126, 4488–4489. (27) Kasson, P. M.; Lindahl, E.; Pande, V. S. PLoS Comput. Biol. 2010, 6, e1000829. (28) Mirjanian, D.; Dickey, A. N.; Hoh, J. H.; Woolf, T. B.; Stevens, M. J. J. Phys. Chem. B 2010, 114, 11061–11068. (29) Sevink, G. J. A.; Zvelindovsky, A. V. Macromolecules 2005, 38, 7502–7513. (30) Lohse, B.; Bolinger, P.-Y.; Stamou, D. J. Am. Chem. Soc. 2008, 130, 14372–14373. (31) Marrink, S. J.; Tieleman, D. P.; Mark, A. E. J. Phys. Chem. B 2000, 104, 12165–12173. (32) Lifshitz, I. L.; Slyozov, V. V. J. Phys. Chem. Solids 1961, 19, 35–50. (33) Chiantia, S.; Ries, J.; Schwille, P. Biochim. Biophys. Acta 2009, 1788, 225–233. (34) Nel, A.; Xia, T.; M€adler, L.; Li, N. Science 2006, 311, 622–627. (35) Jiang, W.; Kim, B. Y. S.; Rutka, J. T.; Chan, W. C. W. Nat. Nanotechnol. 2008, 3, 145–150. (36) Wong-Ekkabut, J.; Baoukina, S.; Triampo, W.; Tang, I.-M.; Tieleman, D. P.; Monticelli, L. Nat. Nanotechnol. 2008, 3, 363–368.

547

dx.doi.org/10.1021/la2033803 |Langmuir 2012, 28, 541–547