A series of molecular dynamics (MD) simulations has been undertaken to investigate the effective interaction between vesicles including PC (phosphatidylcholine) and PE (phosphatidylethanolamine) lipids using the Shinoda–DeVane–Klein coarse-grained force field. No signatures of fusion were detected during MD simulations employing two apposed unilamellar vesicles, each composed of 1512 lipid molecules. Association free energy of the two stable vesicles depends on the lipid composition. The two PC vesicles exhibit a purely repulsive interaction with each other, whereas two PE vesicles show a free energy gain at the contact. A mixed PC/PE (1:1) vesicle shows a higher flexibility having a lower energy barrier on the deformation, which is caused by lipid sorting within each leaflet of the membranes. With a preformed channel or stalk between proximal membranes, PE molecules contribute to stabilize the stalk. The results suggest that the lipid components forming the membrane with a negative spontaneous curvature contribute to stabilize the stalk between two vesicles in contact.
Vesicles made of phospholipids, called “liposomes”, are ubiquitous components of skin moisturizer and other personal-care products and are also utilized for biomedical applications such as drug delivery systems (DDSs) because of their high biocompatibility [1–5]. The stability of vesicles is, however, sometimes an issue in applications, so that several different molecules including water-soluble polymer chains are incorporated to increase the liposome stability [6, 7]. Not surprisingly, lipid composition does significantly affect the stability of liposomes. Understanding the effects of the lipid components on membrane–membrane interactions is likely to be a valuable first step towards the design of stable liposomes on the basis of molecular science. Computer simulations offer a useful way to investigate the effective interaction between membranes from a molecular viewpoint, however, the systems that need to be investigated become too large to be handled with full atomic detail . Especially when the vesicular membranes are used, we have to take into account the role of curvature effects on the interaction. Further, the limitation of sampling efficiency could also be a serious issue because one needs a free energy computation to evaluate the effective interaction between two vesicles .
In this paper, we investigate the effective interaction between vesicles using the recently developed Shinoda–DeVane–Klein (SDK) coarse-grained (CG) force field [8, 10–13]. The SDK-CG model is designed to simulate a membrane with reasonable density, surface tension, molecular partitioning (solvation free energy), as well as distribution functions obtained by all-atom simulations based on the CHARMM force field [14, 15]. In the SDK-CG model, basically, three heavy atoms (and associated hydrogens) are represented by a single CG site. Even though this is a modest CG step, compared with some other CG models [16–18] the computational efficiency is improved by about three orders of magnitude . Thus, the SDK-CG model is useful to investigate the self-assembly processes of lipids and surfactants, vesicle formation, and so on [10, 11, 19–23]. The SDK-CG model has been validated for several different lipid membranes in terms of structure, elastic properties, and tension [12, 24–27]. Previously, we reported that small unilamellar vesicles (SUVs) made of phosphatidylcholine (PC) lipid show high stability against the membrane fusion, with no sign of fusion detected during a several hundred ns-molecular dynamics (MD) simulation even when two vesicles are in contact with each other under an applied external force [9, 12]. This observation suggests that the distance between the centers of mass (COMs) of two vesicles is not a good reaction coordinate to observe membrane fusion. However, taking advantage of this metastability, in this paper, we evaluate the association free energy of two vesicles as a function of distance between them. The effect of lipid composition on the free energy is also investigated in this manner. We also study the stability of a stalk (i.e., a channel linking proximal lipid layers) generated between two apposed vesicles.
We employ the SDK-CG force field to describe the lipid membrane systems . The CG model is designed to reproduce interfacial/surface tension, density, solvation/transfer free energy, as well as distribution functions obtained from all-atomic MD simulations. The CG level and details are given in previous reports [9–12]. We investigate two different lipids: dimyristoyl phosphatidylcholine (DMPC) and dioleoyl phosphatidylethanolamine (DOPE). Vesicles are prepared with 1512 lipids in a manner similar to that used in previous work . First, lipid molecules are randomly disposed in a simulation box with a relatively low degree of hydration, typically about 30 CG water particles per lipid. Typically, a 30 nanosecond (ns) MD trajectory is long enough to observe formation of a continuous aggregated lipid structure in this small system. Even though the structure of this lipid aggregate has not yet converged, nevertheless a random lipid aggregate is achieved at this stage. Then, the so-obtained lipid aggregate is placed in a larger water box to make a dilute lipid solution. Starting from this latter configuration, the lipid aggregate typically transforms into a spherical vesicle within a few hundreds nanoseconds. Three different vesicles have been made by this procedure, composed of DMPC, DOPE, and a 1:1 mixture, respectively.
Molecular dynamics simulations
All MD simulations has been carried out in the NPT ensemble; The Nosé–Hoover thermostat [28, 29] and Parrinello–Rahman barostat  were employed to generate the NPT trajectories. Temperature was controlled at 310 K, and pressure was kept at 0.1 MPa for all MD simulations. The nonbonded interaction, described by LJ functions, was truncated at 1.5 nm, whereas the Coulomb interaction among lipid headgroups was calculated by the PPPM scheme . A time step size of 10 fs was chosen for the MD simulations conducted using the LAMMPS software , in which the SDK-CG model was implemented as the CG/CMM package.
Free energy calculations
Free energy calculations have been carried out by using the thermodynamic integration (TI) scheme . The association energy of two vesicles is defined as the free energy needed to bring them together at a given distance; the reaction coordinate is taken as the distance between the COM of two vesicles, R.
Here the bracket, 〈…〉, denotes the ensemble average over the configurations. kB and T are the Boltzmann constant and temperature, respectively. R0 should be taken at infinitely long distance, though in the actual calculation, R0= 20 nm, where the net interaction between vesicles is zero. The second term in eq. 1 is the Jacobian term, which is subtracted in this representation to compare the free energy profile with that of the vesicle deformation measured separately as shown later (see eq. 3). The constraint MD simulations have been carried out for 100 ns at R′ = 11, 12, 13, 14, 15, 16, 17, 18, 19, and 20 nm, respectively, using a harmonic constraint with the force constant of k = 100 kcal/mol/Å2. This force constant is large enough to fix the two vesicles at the target distance within an error of 0.02 nm. It should be mentioned here that the system containing two small unilamellar vesicles is in a metastable state, and strictly speaking, the calculated free energy along the distance between them is not well defined thermodynamically (in the limit of long time); this is because vesicle fusion potentially takes place, which is not thermodynamically reversible. However, no fusion (or a lipid exchange) was observed during the constrained MD even when two vesicles were in contact with each other (R′ < 16 nm). Taking advantage of the metastability of these vesicles, we use the association free energy defined by eq. 1 to evaluate the vesicle–vesicle interaction. When the distance between vesicles is smaller than the diameter of the vesicle, the calculated energy includes not only net interaction between the vesicles but also a deformation energy due to the contact as shown in Fig. 1b.
The free energy for the vesicle deformation is also evaluated separately. A purely repulsive wall potential on the x–y plane is set at z = z0 in the system. The wall potential is given in a cubic function of z
where K = 100 kcal/mol/Å3, and zi denotes the z-coordinate of the i-th CG particle of a given lipid molecule. This external potential is applied only to lipid molecules, so that water particles can go through the wall without energy penalty. To evaluate the vesicle deformation energy, the COM of a single vesicle is kept at the distance of H from the wall using a harmonic constraint with the force constant of k = 100 kcal/mol/Å2; H = zCOM – z0, where zCOM is the z-coordinate of the vesicle COM. The deformation free energy of vesicle is calculated as a function of H
The constrained MD simulations have been repeated at different distances of H = 5.5, 6.0 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, and 9.5 nm. H0 = 10 nm is safely chosen because no contact of the vesicle with the wall is found in the range of H ≥ 10 nm; ΔG(H) = 0 (H ≥ 10 nm).
Results and discussion
In the present studies, the distance between two vesicles has been kept constant during the CG-MD simulations. However, even with long simulations the vesicles in the present work exhibit high stability against fusion of the two apposed membranes, a finding that differs from the results of previous studies using different CG models [34, 35]. This difference likely arises because membranes made by using the previous CG force field models are known to exhibit spontaneous curvature values that are too negative, which in turn results in a too high propensity for membrane fusion. Specifically, in the present simulations, even though the two vesicles are in contact with each other under an external harmonic potential, no initiation of fusion occurred during a MD trajectory spanning several hundred ns, which effectively corresponds to several μs on a physical time scale (scaled by considering the known enhanced CG-lipid diffusion) . Recently, it was pointed out that the too negative spontaneous curvature problem of previous CG force field models can be improved by changing the lipid model to form membranes with an appropriate spontaneous curvature .This finding might be partly due to the time scale needed for the fusion event. However, it should be pointed out that fusion between “real” biological membranes typically requires fusion machinery or external stimuli [37–40]. The observed metastability of our simulations thus allows changing the distance between COMs of two vesicles, to evaluate the potential mean force (free energy for the association) as a function of the distance between them.
Interestingly, the free energy profile of a pair of SUVs (each with 1512 lipids) is affected by the lipid composition (Fig. 2). DMPC vesicles show a purely repulsive interaction each other (Fig. 2a), although DOPE vesicles show an adhesion with each other at the contact distance of about 16 nm (Fig. 2b). This gives a hint that the DOPE vesicles of this size could be unstable, even though we still did not observe membrane fusion on the sampled time scale. The adhesive nature of PE lipids may explain the fact that a stable vesicle is not attained with only PE lipids experimentally. An adhesive interaction was also observed experimentally for giant unilamellar vesicles . The two contacted vesicles made by a PC/PE mixture showed adhesion only when the vesicle has domain separation. The adhesive interaction is observed through the PE-rich domains, which have negative spontaneous curvature. The vesicles having no such domain separation could not maintain the contact of two vesicles, showing non-adhesive nature . In our MD systems, the DOPE membrane has negative spontaneous curvature and exhibits adhesive interaction between two vesicles.
When two vesicles are kept at a distance of 11 nm, the vesicles are greatly deformed (see Fig. 1b). The free energy penalty at these distances less than the contact distance, i.e., R < 15 nm, may be mostly explained by the energy cost due to the vesicle deformation. To discuss the deformation energy separately, the free energy cost to push the vesicle against a purely repulsive wall is also calculated. In this case, the distance of the COM of the vesicle from the flat repulsive wall is taken as the reaction coordinate of the free energy computation. The deformation free energy in Fig. 2 (dotted lines) is multiplied by 2 in order to make comparison with the two-vesicle system. We also evaluate the net interaction of two vesicles by subtracting the deformation energy from the association free energy of two vesicles. This is also given in Fig. 2. The net interaction between two DMPC vesicles is almost zero, where the repulsive interaction of two vesicles is mostly explained by the deformation energy. In the case of DOPE vesicles, the net interaction is attractive. The deformation energy of DOPE is larger than that of DMPC at the same inter-vesicular distance, because of the larger size of the DOPE vesicle. (DOPE has longer hydrophobic tails by four carbon atoms.) Even though the apposed DOPE membranes show a free energy gain, dehydration of the apposed membranes has not been observed. Figure 3 shows a snapshot of squeezed DOPE vesicles at R = 11 nm, together with water between the vesicles. The water layer is stable enough to give rise to a high free energy barrier for the membrane fusion. Thus, no fusion event is observed during the MD simulations.
The association free energy between two different vesicles made by DMPC and DOPE molecules shows a repulsive interaction (Fig. 2c), although the net interaction seems attractive. In the case of the vesicles made by the 1:1 DMPC/DOPE mixture, a weak adhesion is detected. Indeed, again, the net interaction between the mixed vesicles is attractive. An interesting observation is that the deformation energy of the mixed vesicle is much lower than those of DMPC and DOPE vesicles. Population analysis of lipid species reveals that the higher flexibility of the mixed vesicle originates from the non-uniform distribution of the lipid component within the bilayer.
Figure 4 gives density distribution of lipid components in the DMPC/DOPE mixed vesicles. As in the panel (a), the mixed vesicle has more DMPC and DOPE lipids in the outer and inner leaflets, respectively. This is not a complete segregation of the lipids. As indicated in the method section, we used vesicles transformed spontaneously from a random lipid aggregate in dilute aqueous solution during a CG-MD simulation. The actual numbers of the lipid molecules are 209 DMPC and 300 DOPE lipids in the inner leaflet, respectively, and 547 DMPC and 456 DOPE lipids in the outer leaflet, respectively. To confirm the trend of unequal distribution of lipids, we repeated the CG-MD from 10 different random aggregates and found that the numbers of the lipids are more or less the same within an error of 10 at most. When the vesicles are constrained to a distance of 11 nm, the population of lipids is dramatically changed so as to reduce the elastic energy loss due to the vesicle deformation (Fig. 4b). In the contact area facing to the counter vesicle, the vesicular membrane is flattened. The DOPE lipids, which prefer the negative curvature, are more likely found in the flattened region in the outer leaflet to escape from the positively curved membrane surface. The DMPC lipids in the inner leaflet are populated in the flattened region to prevent the negatively curved surface. No flip-flop motion is detected during the CG-MD simulation so that the domain separation occurs within each leaflet. It should be noted that this is not a separation of the lipid species, but just a population change in a fluidic membrane. The lipids migrate between flattened and curved region, due to lateral diffusion. Again, the flattened outer layer is hydrated as shown in Fig. 3. A similar population change in lipid components is also observed in the distorted vesicle constrained near a flat repulsive wall (Fig. 4c). However, the population change in the flattened membrane region is not so evident, compared with the inter-vesicular compression. This means that the attractive interaction between PE segments in the proximal membranes also help a domain separation to some extent.
As seen above, just by pulling the COMs of two vesicles to contact with each other, we could not observe even an initiation of fusion, which is believed to start as stalk formation between the proximal membranes [42, 43]. Here, we have investigated the stability of a stalk by a CG-MD simulation started from two vesicles connected by a stalk-like channel. The initial channel, as shown in the left top panel of Fig. 5, has a bilayer thickness, which connects the outer leaflets of two vesicles. The stability of the channel is changed largely depending on the lipid components. The channel in the DMPC vesicles gradually thins and is eventually eliminated due to thermal fluctuations under the repulsive nature of the vesicle–vesicle interaction. This happens within 100 ns MD, though the lifetime of the channel may depend on the initial channel structure. On the contrary, the channel between the DOPE vesicles thickens gradually, keeping the bilayer thickness in the other direction, and grew circular on the vesicle surface to increase the contact between the two vesicles (see Fig. 5 lower middle panel). As a result, an inverse (elongated) micelle-like structure is formed between the two vesicles. It should be noted that the inverse structure is preferable for DOPE molecules under dehydration. This transformation has happened within 100 ns. Another 100 ns CG-MD shows a convergence of the morphology of the connected vesicles. This is totally different from an ordinary stalk. The inner leaflets of different vesicles have no contact with each other, though the lipids in the outer leaflets diffuse across the vesicles. It should be pointed out that the structure observed here is similar to the hemifusion model proposed to explain the strong adhesive behavior between vesicles containing a domain of lipid having a negative spontaneous curvature . In the case of the 1:1 DMPC/DOPE vesicles, the channel stays stable with the thickness of the bilayer. Extending the CG-MD run up to 700 ns confirms no change in the morphology. A population analysis of lipids reveals that DOPE molecules are more likely found around the channel to support the negatively curved surface of the proximal membrane. The channel is the only portion of the outer membrane having the negative curvature. This is clearly seen in the density map given in the top right panel in Fig. 5. This result suggests that PE lipids contribute to increase the stability of the stalk structure. Thus, mixed membranes adjust their inner structure in the response to either topological or morphological changes, a process that contributes to reduce the free energy cost required for any imposed structural changes of the lipid membranes.
We have carried out CG-MD simulations to investigate the stability of vesicles composed of PC and PE lipids. We investigated only small vesicles made from 1512 lipid molecules, and find that two apposed vesicles are stable against membrane fusion. This observation holds even when the two vesicles are contacted to each other under an applied external harmonic potential. A water layer stays between the vesicles even when the pulling force is large enough to deform the vesicles. This may suggest that headgroup dehydration is a key step to observe membrane fusion. The association free energy of the two vesicles computed by thermodynamic integration depends on the lipid composition of the vesicles. Two DMPC vesicles are purely repulsive to each other, though DOPE vesicles show a free energy gain at the contact distance. A mixed vesicle of DMPC and DOPE (1:1) shows a higher flexibility having a lower energy barrier on compression, which is caused by a two-dimensional population change between the two lipid components in each leaflet of the membrane. With a generated channel between proximal membranes, DOPE molecules, which form a membrane with negative spontaneous curvature, contributed to stabilize the channel. The present results suggest that the lipid components forming the membrane with a negative spontaneous curvature contribute to stabilize the stalk between two vesicles in contact.
A collection of invited papers based on presentations at the 33rd International Conference on Solution Chemistry (ICSC-33), Kyoto, Japan, 7–12 July 2013.
This research was supported by JSPS KAKENHI Grant Number 23350014, Next Generation Supercomputing Project, Nanoscience Program, TCCI/CMSI in the Strategic Programs for Innovative Research, MEXT, Japan, and HPCI Systems Research Projects (Project ID hp120093). The numerical calculations were mainly carried out on the T2K Open Supercomputer at the Center for Computational Sciences, University of Tsukuba, the TSUBAME2.0 supercomputer in the Tokyo Institute of Technology, and the supercomputer of ACCMS, Kyoto University. MLK thanks NSF Grant 1212416 – Building Computational Models to Probe Membrane Fusion – for partial support as well as the Temple High-Performance GPU/CPU System that was purchased in part with funds from NSF Grant MRI-R2 095885.
 T. M. Allen, P. R. Cullis. Adv. Drug. Deliv. Rev. 65, 36 (2013).Search in Google Scholar
 T. Lian, R. Ho. J. Pharm. Sci. 90, 667 (2001).Search in Google Scholar
 A. Samad, Y. Sultana, M. Aqil. Curr. Drug. Deliv. 4, 297 (2007).Search in Google Scholar
 J. de Leeuw, H. C. de Vijlder, P. Bjerring, H. A. M. Neumann. J. Eur. Acad. Dermatol. Venereol. 23, 505 (2009).Search in Google Scholar
 D. Lasic, D. Needham. Chem. Rev. 95, 2601 (1995).Search in Google Scholar
 W. Shinoda, R. DeVane, M. L. Klein. Curr. Opin. Struct. Biol. 22, 175 (2012).Search in Google Scholar
 W. Shinoda, R. Devane, M. L. Klein. Mol. Simulat. 33, 27 (2007).Search in Google Scholar
 A. MacKerell Jr., D. Bashford, M. Bellott, R. L. Dunbrack Jr., J. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher III, B. Roux, M. Schlenkrich, J. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiorkiewicz-Kuczera, D. Yin, M. Karplus. J. Phys. Chem. B 102, 3586 (1998).10.1021/jp973084fSearch in Google Scholar PubMed
 Z. Wu, Q. Cui, A. Yethiraj. J. Chem. Theory Comput. 7, 3793 (2011).Search in Google Scholar
 M. Venturoli, M. M. Sperotto, M. Kranenburg, B. Smit. Phys. Rep. 437, 1 (2006).Search in Google Scholar
 R. DeVane, W. Shinoda, P. B. Moore, M. L. Klein. J. Chem. Theory Comput. 5, 2115 (2009).Search in Google Scholar
 B. G. Levine, D. N. LeBard, R. DeVane, W. Shinoda, A. Kohlmeyer, M. L. Klein. J. Chem. Theory Comput. 7, 4135 (2011).Search in Google Scholar
 T. Nakamura, W. Shinoda, T. Ikeshoji. J. Chem. Phys. 135, 094106 (2011).Search in Google Scholar
 T. Nakamura, W. Shinoda. J. Chem. Phys. 138, 124903 (2013).Search in Google Scholar
 K. Shuhei, T. Nakamura, S. O. Nielsen, W. Shinoda. J. Chem. Phys. 139, 034108 (2013).Search in Google Scholar
 S. Nosé, M. L. Klein. Mol. Phys. 50, 1055 (1983).Search in Google Scholar
 M. Parrinello, A. Rahman. J. Appl. Phys. 52, 7182 (1981).Search in Google Scholar
 S. Marrin, A. Mark. J. Am. Chem. Soc. 125, 11144 (2003).Search in Google Scholar
 M. Stevens, J. Hoh, T. Woolf. Phys. Rev. Lett. 91, 188102 (2003).Search in Google Scholar
 S. J. Marrink, M. Fuhrmans, H. J. Risselada, X. Periole. In Coarse-Graining of Condensed Phase and Biomolecular Systems, G. A. Voth (Ed.), pp. 5–19, CRC Press, Boca Raton (2008).Search in Google Scholar
 J. Rizo, C. Rosenmund. Nat. Struct. Mol. Biol. 15, 665 (2008).Search in Google Scholar
©2014 IUPAC & De Gruyter Berlin Boston