Coarse-grained molecular dynamics simulation of binary charged lipid membranes: Phase separation and morphological dynamics

Coarse-grained molecular dynamics simulation
of binary charged lipid membranes:
Phase separation and morphological dynamics

Hiroaki Ito Department of Mechanical Engineering, Graduate School of Engineering, Osaka University, Osaka 565-0871, Japan    Yuji Higuchi Institute for Materials Research, Tohoku University, Miyagi 980-8577, Japan    Naofumi Shimokawa School of Materials Science, Japan Advanced Institute of Science and Technology, Ishikawa 923-1292, Japan
July 15, 2019

Biomembranes, which are mainly composed of neutral and charged lipids, exhibit a large variety of functional structures and dynamics. Here, we report a coarse-grained molecular dynamics (MD) simulation of the phase separation and morphological dynamics in charged lipid bilayer vesicles. The screened long-range electrostatic repulsion among charged head groups delays or inhibits the lateral phase separation in charged vesicles compared with neutral vesicles, suggesting the transition of the phase-separation mechanism from spinodal decomposition to nucleation or homogeneous dispersion. Moreover, the electrostatic repulsion causes morphological changes, such as pore formation, and further transformations into disk, string, and bicelle structures, which are spatiotemporally coupled to the lateral segregation of charged lipids. Based on our coarse-grained MD simulation, we propose a plausible mechanism of pore formation at the molecular level. The pore formation in a charged-lipid-rich domain is initiated by the prior disturbance of the local molecular orientation in the domain.

I Introduction

Biomembranes, which are formed by self-assembly of various amphiphilic molecules, exhibit a large variety of functional structures and dynamics because these membranes are small, soft, and heterogeneous. There is dynamic compositional heterogeneity in biomembranes, namely, the transient formation of small domains enriched in saturated lipids and cholesterol, known as lipid rafts Eggeling et al. (2009); Toulmay and Prinz (2013). Raft domains may play important roles in signal transduction and membrane trafficking Simons and Ikonen (1997); Simons and Sampaio (2011); Suzuki et al. (2012). To study the physicochemical mechanism of this lateral membrane heterogeneity, phase separation in multi-component lipid membranes, especially in giant unilamellar vesicles (GUVs) consisting of mixtures of saturated and unsaturated lipids, has attracted a great deal of attention as a model in vitro.

Typically, in ternary lipid mixtures composed of saturated lipids, unsaturated lipids, and cholesterol, liquid-ordered () and liquid-disordered () phases coexist Veatch and Keller (2002); Baumgart et al. (2003); Komura and Andelman (2014). In particular, the equilibrium conditions for stable phase separation Veatch and Keller (2003); Heberle and Feigenson (2011) and the domain-growth dynamics after the sudden transition into coexistence conditions Saeki et al. (2006); Yanagisawa et al. (2007); Stanich et al. (2013) have been extensively investigated, confirming both the equilibrium and dynamic features of phase separation in electrically neutral systems. Recently, there has been increasing interest in the effects of long-range electrostatic interactions on the phase separation Vequi-Suplicy et al. (2010); Shimokawa et al. (2010); Blosser et al. (2013); Pataraia et al. (2014); Himeno et al. (2014) because plasma membranes van Deenen et al. (1982); Yeung et al. (2008) and some organelles, such as lysosomes Dowhan (1997) and mitochondria Schlame (2008), are enriched with charged lipids. In experiments with systems containing charged lipids, charged unsaturated lipids tend to suppress the phase separation compared with neutral lipid systems Vequi-Suplicy et al. (2010); Shimokawa et al. (2010); Blosser et al. (2013); Pataraia et al. (2014); Himeno et al. (2014), whereas charged saturated lipids promote phase separation Himeno et al. (2014). Several theoretical models May et al. (2002); Shimokawa et al. (2011); Okamoto et al. (2016); Shimokawa et al. (2016) have also been proposed to explain the electrostatic effects on the lateral phase separation in charged lipid membranes.

Morphological changes in lipid bilayer membranes are important dynamic phenomena in biomembranes. Examples include the budding in endo- and exocytosis McMahon and Gallop (2005), the fission-fusion sequence of vesicular transport Beckers et al. (1989), and pore formation in autophagy Juhasz and Neufeld (2006). Inspired by the dynamic changes in biomembrane morphology, studies have also investigated the morphological changes in GUVs induced by external stimuli, such as osmotic pressure Boroske et al. (1981), addition of surfactant molecules Hamada et al. (2009), and a synthetic photosensitive amphiphile Hamada et al. (2010). Some ions or charged biomolecules can also trigger dynamic morphological changes through adsorption onto the lipid monolayers or bilayers via the electrostatic interactions. For example, locally injected HCl drives cristae-like deformation Khalifat et al. (2008), BAR (Bin/amphiphysin/Rvs) domain family proteins induce the formation of tubular structures on the membranes Saarikangas et al. (2009), and actin cytoskeletons with myosin motor proteins cause concave deformations of oppositely charged lipid membranes Ito et al. (2015); Nishigami et al. (2016).

In biomembranes, phase separation and morphological changes in the membrane are closely coupled. Recently, raft domains have been reported to play a critical role in membrane dynamics, such as autophagy and endocytosis, in response to chemical signals Rajendran and Simons (2005); Dall’Armi et al. (2013). Using GUVs composed of neutral lipids, coupling between phase separation and morphological dynamics has also been observed. For example, the adhesion between GUVs consisting of inverse cone- and cylinder-shaped lipids Sakuma et al. (2008), the pore formation in membranes consisting of cone- and cylinder-shaped lipids Sakuma et al. (2010), and the coupling between the phase separation and deformed membranes by osmotic pressure Yanagisawa et al. (2008); Hamada et al. (2007) have been reported. Several studies have focused on the coupling between the phase separation and the morphology in charged lipid membranes. For instance, a positively charged protein, cytochrome c, induces the collapse of domains enriched with cardiolipin (CL), which is a negatively charged unsaturated lipid abundant in the inner membrane of mitochondria Beales et al. (2011). Even for more physiologically relevant phospholipids with ions, the domains enriched with the unsaturated negatively charged lipid (dioleoylphosphatidylserine; DOPS) bud toward the interior of the GUV upon the addition of calcium ions Shimokawa et al. (2010), and those enriched with the saturated charged lipids (dioleoylphosphatidylglycerol; DOPG) undergo pore formation in the membrane by electrostatic repulsion Himeno et al. (2015).

Although phase separation and the morphologies of charged lipid membranes have been extensively investigated, there are still several questions about the mechanism of charge-specific dynamics. For example, although the domain-coarsening dynamics are frequently examined in neutral membranes Saeki et al. (2006); Yanagisawa et al. (2007); Stanich et al. (2013), these experiments have not been performed in charged lipid membranes due to the difficulty in preparing experimental systems in which the charge and ion conditions are well-controlled. It is difficult to obtain reliable statistics because the established vesicle-preparation method of electroformation is not applicable to charged systems Blosser et al. (2013). To reveal the phase-separation dynamics, morphological changes, and the coupling between them in charged lipid membranes, we use a coarse-grained molecular dynamics (MD) simulation Markvoort et al. (2009); Zheng et al. (2010) of binary charged and neutral lipid mixtures. To capture such mesoscopic and long-time features of vesicle dynamics, highly coarse-grained model is necessary Venturoli et al. (2006); Noguchi (2009), as we successfully reproduced the experimental morphologies of charged lipid vesicles in our previous work using the present simulation model Himeno et al. (2015). The coarse-grained MD simulation is suitable for this purpose because we can specify various conditions, such as lipid composition and salt concentration, we can easily investigate the effects of pure electrostatic interactions on a membrane composed of a number of dynamically moving molecules, and we can observe the underlying mechanism at an adequate molecular level. We investigate the dynamic growth of the charged domain and the coupling between the phase separation and the morphological changes including topological changes. The details of the coarse-grained MD simulation are introduced in Sec. II. The phase separation and the domain growth are investigated in Secs. III.1 and  III.2. The morphological changes induced by electrostatic interactions are presented in Sec. III.3. The implications of our work and comparisons with other studies are discussed in Sec. IV.

Ii Methods

We extended the coarse-grained model described by Cooke et al. Cooke et al. (2005) to calculate the behavior of charged lipid membranes. This highly coarse-grained model has been used to simulate various mesoscopic behaviors Cooke et al. (2005); Guigas and Weiss (2006), where the physical properties of the same orders of magnitude with more detailed models have been obtained Wang and Deserno (2016). A lipid molecule is represented by one hydrophilic head bead followed by two hydrophobic tail beads. The excluded-volume interaction between two beads separated by a distance is


where . is the unit of energy. We chose and to form a stable bilayer, where Å is the unit of length corresponding to the cross-sectional diameter of a single lipid molecule. The potentials for the stretching and bending of bonds between connected beads are expressed as




where and are the bonding strength of connected beads and the bending stiffness of a lipid molecule, respectively. is the angle between adjacent bond vectors. The hydrophobic attractive interaction among hydrophobic beads is expressed as


where is the cut-off length for the attractive potential. The lipid membrane is in the gel phase when is large, whereas it is in the liquid phase when is small. Such gel and fluid phases observed in the binary lipid bilayer membranes are directly relevant to the liquid-ordered () and liquid-disordered () phases, respectively, which are experimentally observed in ternary lipid bilayer membranes composed of unsaturated lipids, saturated lipids, and cholesterols, as both lateral phase separations in the simulation and experiments are recognized by circular domains. Thus, it is possible to determine whether the lipid species is in the gel or liquid phase by controlling in a neutral lipid system Cooke et al. (2005). To represent the charged lipids, we considered the electrostatic interaction among charged head groups in addition to the other interactions. The electrostatic repulsion is described as the Debye-Hückel potential


where is the Bjerrum length, and are the valencies of the interacting charged head groups, and is the Debye screening length which is related to the bulk salt concentration (, , , and are the dielectric constant of the solution, Boltzmann constant, absolute temperature, and elementary charge, respectively). To capture the behaviors caused by the screened electrostatic interaction qualitatively, is set to , , or , varying by an order of magnitude around the physiological salt concentration of . Because the typical negatively charged head groups, such as phosphatidylglycerol (PG), phosphatidylserine (PS), phosphatidylinositol (PI), and phosphatidyl acid (PA), have a monovalent charge, we set . We did not set a cutoff for this screened electrostatic interaction.

Each bead position obeys the stochastic dynamics described by the Langevin equation


where and are the mass and drag coefficients, respectively. The force is calculated from the derivatives of the interaction potentials Eqs.(1)–(5). The constant is chosen as the unit for the timescale, and the time increment is set at . The Brownian force satisfies the fluctuation-dissipation theorem


In this study, we calculated bilayer vesicles composed of two lipid species to represent the ternary lipid vesicles with gel and liquid domains observed in experiments Cooke et al. (2005). The calculated lipid composition is a binary mixture of A-lipids with or without a monovalent electric charge at the head beads, corresponding to charged or neutral lipids, respectively, and neutral B-lipids. Binary neutral membranes consisting of neutral A- and B-lipids are used as the reference condition for the binary charged membrane consisting of charged A-lipids and neutral B-lipids. The initial state of the calculation is a spherical lipid bilayer consisting of 5000 lipid molecules in total, where A- and B-lipids are mixed homogeneously at a specified ratio, A:B. The total calculation time is set to , during which the phase-separation or morphological dynamics adequately relaxes.

Iii Results

iii.1 Lateral phase separation

Figure 1: (a) Typical snapshots of the lateral phase separation of a neutral vesicle consisting of 2500 neutral A-lipids [head, red (dark); tail, green] and 2500 neutral B-lipids [head, yellow (light); tail, cyan]. After rapid relaxation on the order of , the binary vesicle reaches a quasi-equilibrium state. Attractive interaction parameters are set to , , and . (b) Time evolution of the attractive potential per lipid molecule during the lateral phase separation for under various attractive interactions between A-lipids of , , , , and . Calculations were performed five times for each value to ensure reproducibility.
Figure 2: (a) Mean-squared diffusion length of a lipid in plotted versus attractive interaction parameter . (b) Nematic order parameter plotted versus . The values are averaged spatially over a vesicle and temporally over for each in both (a) and (b). Error bars represent the standard deviations. Other interaction parameters are set to , and .

First, we investigated the dynamic features of lateral phase separation in neutral binary lipid vesicles as the reference behavior for that of charged lipid vesicles. Because both A- and B-lipids are neutral, we set the valencies of their head groups as . In other words, force is obtained from Eqs.(1)–(4) without the electrostatic potential in Eq.(5). Figure 1(a) shows typical snapshots of lateral phase separation in the binary mixture of neutral A-lipids (red or dark) and B-lipids (yellow or light) with A:B = 2500:2500, where the attractive interaction parameters were set to , , and . Immediately after the start of the simulation, the vesicle was in the laterally homogeneous state (). The phase separation was gradually induced by the attractive interaction being stronger among the same species than that among different species on a timescale of . Figure 1(b) shows the attractive potential, , per lipid molecule during the phase separation for attraction parameters , , , , and , where the bilayer membrane phase ranged from fluid to gel Cooke et al. (2005). We quantitatively confirmed the transition between liquid and gel phases induced by the attractive interaction parameter at , based on the calculations of the lipid diffusivity and nematic order parameter in the vesicle. Figure 2(a) shows the mean-squared diffusion length of a lipid in plotted versus . Although a smaller is accompanied by a larger error resulting from spike noises caused by frequent dropouts of lipids, the diffusion length decreased according to the increase in , and reached a plateau for . We also calculated the nematic order parameter by , where is the unit bond vector of tail beads of the th lipid and represents the average over all the combinations of lipid pairs [Fig. 2(b)]. The nematic order parameter also shows a clear plateau at , indicating the most ordered state with the lowest diffusivity in this parameter range. These two features in Figs. 2(a) and 2(b) characterize the transitional behavior of lipid bilayer membranes between gel and liquid phases in our simulation. Despite the phase transition around , even though the values of were slightly different according to in Fig. 1(b), the vesicles exhibited the same qualitative feature in their dynamics, namely, a decrease in potential on a timescale of . After this rapid relaxation into phase-separated states, the vesicles reached quasi-equilibrium states, identified by the flattening of the temporal changes in . Complete phase separation with two macro domains is reached stochastically on much longer timescales than the practical calculation time in this study because of the much slower diffusion of large domains compared with single molecule. Because the other kinds of energies , , and did not exhibit any clear changes during the phase separation, the change in attractive interaction dominated the lateral phase separation in neutral lipid vesicles in both the gel and liquid states.

Figure 3: (a) Typical snapshots of the lateral phase separation of a charged vesicle consisting of 2500 charged A-lipids [head, red (dark); tail, green] and 2500 neutral B-lipids [head, yellow (light); tail, cyan] at salt concentrations and . Attractive interaction parameters are set to , , and for both salt conditions. (b) Attractive and (c) electrostatic repulsive potentials and , respectively, for under various attractive interactions between A-lipids of , , , , and at . (d) Attractive and (e) electrostatic repulsive potentials at . Calculations were performed three times for each and value to ensure reproducibility.

Next, we considered a negatively charged vesicle consisting of negatively charged A-lipids and neutral B-lipids to examine the effect of long-range electrostatic repulsion on the lateral phase separation. The valencies of the head groups of A- and B-lipids were set as and , respectively. Figure 3(a) shows typical snapshots of the lateral phase separation of vesicles composed of 2500 charged A-lipids and 2500 neutral B-lipids under different salt concentrations. All the parameters except for the valency of the A-lipids were the same as in the calculation for the neutral vesicles (Fig. 1).

Figure 4: (a) Typical snapshots of the lateral phase separation of neutral (upper) and charged (at , lower) vesicles for , , , , and . Other conditions are set to , A:B = 2500:2500. (b) Corresponding attractive potentials . Calculations were performed three times for each charge condition and value to ensure reproducibility.

For charged lipids at a higher salt concentration of the vesicles exhibited dynamics similar to those of the neutral vesicles, whereas at a lower salt concentration of they exhibited clearly slower phase-separation dynamics. The attractive potential of the charged lipid vesicles at (Fig. 3(b)) was also quantitatively the same as that of the neutral lipids (Fig. 1(b)), although the electrostatic repulsive potential, , increased with the phase-separation process (Fig. 3(c)). Thus, the electrostatic interaction among charged lipids with was negligible at , and the increase in was two orders of magnitude lower than the decrease in . However, at , the phase-separation dynamics were affected by the less-screened electrostatic repulsion as the decrease in became slower and increased by an order of magnitude compared with that at (Fig. 3(d) and  3(e)).

We compared the static phase-separated states of neutral and charged vesicles consisting of 2500 neutral or charged A-lipids and 2500 neutral B-lipids. We focused on the dependence of the static states because the attractive parameter governs the degree of mixing, which provides quantitative information about susceptibility of the phase separation in neutral and charged vesicles. As shown in Fig. 4(a), there was a large difference in the the critical attractive parameter, , for the lateral phase separation. The corresponding changes in the attractive interaction potential, , shown in Fig. 4(b) clearly indicate the transition point of , which coincides with the constant interaction potential over . The plateau in occurred at for neutral vesicles, but at for charged vesicles at . This inhibition of the lateral phase separation by the electrostatic repulsion among charged lipids qualitatively agrees with previous experiments Vequi-Suplicy et al. (2010); Shimokawa et al. (2010); Blosser et al. (2013); Pataraia et al. (2014); Himeno et al. (2014) and theoretical calculations based on the free energy of charged membranes Shimokawa et al. (2010); May et al. (2002).

iii.2 Coarsening dynamics

Figure 5: Density correlation functions between lipid species. The intersection between the (A-A) correlation curve and (A-B) correlation curve is the typical domain size of A-lipid-rich domains. Typical conditions, where , , A:B=1500:3500, and , are used as an example.
Figure 6: (a) Time evolution of domain size for neutral vesicles consisting of neutral A- and B-lipids (black filled circles) and charged vesicles consisting of charged A-lipids and neutral B-lipids at (purple filled squares) and (red open squares). Lipid composition of A:B = 2500:2500. Error bars represent the standard deviation at each time. Calculations were performed five times for neutral and charged lipids at , and three times for charged lipids at . Slopes for and in the double logarithmic axis are shown as visual guides. (b) Lipid composition of A:B = 1500:3500 in the same manner in (a). Calculations were performed twenty times for neutral, and ten times for charged lipids at and . Slopes for , , and are shown as visual guides.

The dynamics of phase separation in lipid vesicles has been extensively studied experimentally and theoretically in terms of a unified view of scaling behaviors of the domain coarsening, where the domain size is proportional to  Stanich et al. (2013); Laradji and Kumar (2005) and is the domain-growth exponent. Thus, we analyzed for neutral and charged lipid vesicles. To define the typical domain size, the density correlation function averaged over the vesicle was calculated. Figure 5 shows the spatial distance at which the density autocorrelation function of A-lipids (A-A) intersects with the density cross-correlation function between A- and B-lipids (A-B); this distance is defined here as the typical size of A-lipid-rich domains. Figure 6 shows the different coarsening dynamics in the neutral lipid vesicles, charged lipid vesicles for , and charged lipid vesicles for at the lipid compositions of A:B = 2500:2500 (Fig. 6(a)) and 1500:3500 (Fig. 6(b)). At A:B = 2500:2500, the domain-growth exponents of neutral A-lipid-rich domains and charged A-lipid-rich domains at were both (Fig. 6(a)), which is generally seen in the coarsening process caused by domain-area deformation driven by line tension Saeki et al. (2006); Yanagisawa et al. (2007); Stanich et al. (2013). Under these conditions, the observed phase separation occurred by spinodal decomposition with percolated domain. This bicontinuous structure with a clear domain edge was deformed over time by the line tension between the A- and B-lipid-rich phases, which originated from the different attractive interactions between the same and different species. In contrast, the coarsening of the charged A-lipid-rich domains at showed slower dynamics and , which is generally seen during evaporation condensation (Ostwald ripening), where the growth exponent of is referred to as the Lifshitz–Slyozov–Wagner model Lifshitz and Slyozov (1961); Wagner (1961), or during diffusion-limited collision and coalescence of small liquid domains Laradji and Kumar (2005). Figure 3(a) at exemplifies the different slow coarsening dynamics in the binary vesicles in the presence of relatively strong long-range repulsion. The domain edge was blurred in this case, and a number of small domains randomly came together rather than a single percolated domain undergoing tension-driven deformation .

In addition, we examined the lipid composition of A:B = 1500:3500 (Fig. 6(b)), for which the area fraction of charged A-lipids was almost , because the coarsening exponent depends on the area fraction of the phase-separated domains Stanich et al. (2013). The coarsening dynamics of the neutral lipid vesicles and charged lipid vesicles at were again similar. After the initial domain formation with a domain size of more than 2 or 3 lipids, the domain-growth exponent was as in the case of A:B = 2500:2500; however, decreased to or slightly less, probably due to the small number of A-lipids and A-lipid-rich domains. The small number of lipids formed small circular domains rather than a percolated domain, resulting in the collision and coalescence of the liquid domains with an exponent of . The coarsening dynamics of charged lipid vesicles at were also an order of magnitude slower than those under neutral or well-screened charged conditions. Interestingly, the vesicles remained for a time during the initial homogeneous state and then suddenly started coarsening, with an ill-defined exponent of or less.

iii.3 Morphological changes

Figure 7: (a) Typical snapshots of the morphological changes of a binary charged vesicle composed of 2500 charged A-lipids (head, red (dark); tail, green) and 2500 neutral B-lipids (head, yellow (light); tail, cyan) at . , , and . The upper and lower image series respectively show surface and cross-sectional images of the same vesicle. Typically, the morphological changes occur an order of magnitude slower () than the lateral phase separation (). Time evolution of the (b) attractive and (c) electrostatic repulsive potentials and per lipid molecule, respectively. The calculation was performed once and showed plausible results (see Fig. 8).

In this section, we turn our focus to the morphological dynamics of the charged lipid vesicles. In the previous section, we calculated the phase-separation dynamics on a timescale of . However, when we continued the calculation until , we observed the morphological changes in the charged lipid vesicles depending on the salt concentration . In contrast, we did not observe morphological changes within for neutral lipid vesicles (data not shown). Figure 7(a) shows typical snapshots of the morphological changes, including the topological change from spherical to disk-shaped, which were observed during the longer time calculation at . The lipid composition was set to A:B = 2500:2500, and the interaction parameters were set to , , and . After the relaxation of the lateral phase separation for (second image in Fig. 7(a)), a pore was formed at a charged-lipid-rich domain (red or dark) in the membrane at around . The pore diameter increased with time and finally the vesicle transformed into a disk shape. The charged lipids were segregated at the edge of the disk, whereas the neutral lipids formed a planar bilayer membrane surrounded by the charged lipids (last image in Fig. 7(a)). During the pore formation, the almost relaxed attractive potential began to increase (Fig. 7(b)), whereas the increased electrostatic repulsive potential started to decrease (Fig. 7(c)). After the onset of pore formation at , these potentials changed as the vesicle morphology relaxed to the flat disk shape in . Although the potentials slightly changed even after , due to the gradual dropouts of lipid molecules from the structure, we confirmed that the disk shape was kept up to . We thus regarded the morphologies at as the static states of the pre-assembled vesicles. Note that such states observed in experiments have been treated as practically static because of the much slower timescale of complete “melting” of the structure.

Figure 8: (a) Typical morphologies for various salt conditions and attractive interactions between A-lipids. is from to and is from to . and were fixed as and , respectively. (b) Phase diagram of the charged vesicle morphologies. We defined four types of morphologies; spherical vesicle (circle), disk (square), string (triangle), and bicelle (cross). The dashed lines are visual guides indicating the phase boundaries.

Finally, we determined the phase diagram of the vesicle morphology for various salt concentrations and attractive interaction parameters between A-lipids (Fig. 8). The salt concentration, , was varied from to , ranging from almost neutral conditions with a corresponding Debye screening length of (Å) to highly repulsive conditions with (Å). The attraction parameter for A-lipids, , was varied from (liquid state) to (gel state). In this wide parameter region, Fig. 8(a) shows the typical morphologies of charged lipid vesicles coupled to the lateral phase separation observed at , when the vesicles reach static states as shown in Figs. 7(b) and  7(c). At , the lipid vesicles remained spherical and macro-phase separation was observed for all the values of . At , the vesicles exhibited various morphologies depending on . For , the initially spherical vesicles became string-shaped. For , , and , the vesicles became disk-shaped, and for , they remained spherical. At , where the electrostatic repulsion was substantially stronger than the attractive interaction among lipid molecules, the self-assembled structures were no longer maintained. From to or , the vesicles broke into small aggregations, called bicelles, whereas at the vesicles adopted a spiky string morphology. When the transitions to disk, string, or bicelle occurred, the charged lipids assembled at the edge of these static structures. Figure 8(b) shows a phase diagram summarizing these static shapes.

Iv Discussion

In this study, we investigated the effects of screened but still long-range electrostatic repulsion on the dynamic behaviors of lipid bilayer vesicles, especially the lateral-phase-separation dynamics and morphological changes. After a temperature quench below the miscibility transition temperature, the phase-separation dynamics of a homogeneous neutral lipid membrane are nucleation growth from a metastable state or spinodal decomposition from an unstable state. During the nucleation process, there is an effective energy barrier against the localization of a lipid species, owing to factors such as the mixing entropy of lipids and the line tension of domains. In contrast, during spinodal decomposition, a binary lipid mixture immediately segregates into two coexisting phases. The occurrence of these two instability processes generally depends on the level of the parameter quench at the onset of the phase separation and the lipid composition of the vesicle, namely, whether the system state is located below the spinodal line or between the binodal and spinodal lines in the two-parameter field Stanich et al. (2013); Lipowsky and Dimova (2003). For example, Veatch and Keller reported the instability of stripe-shaped domains with only neutral lipids, suggesting spinodal decomposition in a vesicle of (dioleoylphosphatidylcholine; DOPC):(dipalmitoylphosphatidylcholine; DPPC) (1:1) and 35% cholesterol caused by a temperature quench near the critical miscibility temperature Veatch and Keller (2003). After the onset of the phase separation by nucleation or spinodal decomposition, the emergent domains grow further by evaporation condensation (Ostwalt ripening), collision and coalescence by diffusion, and domain deformation by line tension, depending on the area fraction of the domain-forming lipids Stanich et al. (2013). These domain growth processes at this stage are called “coarsening”.

Figure 9: Example of the molecular mechanism of pore formation. The vesicle is composed of 2500 charged A-lipids (head, red (dark); tail, green) and 2500 neutral B-lipids (head, yellow (light); tail, cyan) at . and . Sequential snapshots in the middle show the surface and cross-sectional images taken at to . Dotted squares contain magnified images of the pore-forming site at (before), (onset), and (after), for surface (upper) and cross-sectional (lower) views of the bilayer membrane. Embedded lipids and the surrounding lipids are displayed with spheres of the original size and small spheres, respectively.

The present numerical calculations revealed that the attraction parameters, the lipid composition, and the electrostatic long-range repulsion play important roles in the coarsening process and the membrane instability. Although it is difficult to observe the exact onset of the phase separation owing to the short length and fast timescale, our MD simulation can shed light on this phenomenon. In the coarsening dynamics of charged lipid membranes at for the lipid composition A:B = 2500:2500, the lateral structure of the membrane remained homogeneous for a certain duration, whereas the neutral lipid membranes or charged lipid membranes at immediately began the phase separation (Fig. 6(b)), even though the parameter conditions were the same. This result strongly indicates that the dynamics of the charged lipid membranes at the onset of the phase separation are not spinodal decomposition, but instead nucleation growth starting from a metastable state with an energy barrier, probably originating from the electrostatic repulsive force among charged lipids, that needs an overcoming time . During late phase separation, the domain growth in charged membranes was much smaller than that in neutral membranes (Fig. 6), and growth exponent values were also different from those in neutral membranes. Pataraia et al. previously observed transient small circular domains or bicontinuous domains even after 20 min from vesicle preparation of ternary charged vesicles composed of negatively charged dioleoylphosphatidylglycerol (DOPG), egg sphingomyelin (eSM), and cholesterol, indicating the slower coarsening dynamics. This finding agrees well with our numerical calculations and implies the underlying general physicochemical effects on the phase-separation process of electrostatic interaction among charged lipids. However, the details of the coarsening dynamics, such as growth exponent of charged lipid membranes, remains less well understood. In this study, we observed an ill-defined exponent for A:B = 1500:3500 [Fig. 6(b)]. Two-dimensional charged membranes can be exempted from electroneutrality because their counterions can escape into the three-dimensional bulk, and so the coarsening dynamics in charged membranes might be a unique feature of a two-dimensional electrolyte solution. Therefore, examining for charged lipid membranes in other experimental or theoretical works would be an interesting step toward understanding the dynamic heterogeneity in charged membrane systems.

Li et al. recently reported the first systematic theoretical study of the equilibrium conformational changes caused by the electrostatic interactions among surface lipids by calculating the minimum bending and Coulombic energies Li et al. (2015). Although they assumed a fixed uniform surface charge distribution, our model can treat a more plausible situation, in which the discrete charged and neutral lipid mixtures move laterally in the membrane. Thus, we found the characteristic topological change caused by pore formation in a charged membrane (see Fig. 8). In addition, our MD simulation visualized details of the behavior at a molecular level. Figure 9 shows an example of the dynamics of lipid molecules around the pore-forming site within the charged-lipid-rich domain (red or dark) during the topological change. Interestingly, even before the onset of the pore formation at , the orientation of several charged lipids became perpendicular to that of the surrounding lipids at . The charged head group of the perpendicular lipid was embedded between the two leaflets of the bilayer membrane, probably because of the strong electrostatic repulsion among their charged head groups. This could trigger further instability in the membrane structure by interfering with the hydrophobic ordering attraction of the surrounding lipid chains, resulting in the pore formation.

It should be noted that in the present calculation, the vesicle diameter was on the order of due to the limitation imposed by computational cost. Previous studies have also performed similarly sized coarse-grained simulations to reproduce the dynamics of neutral lipid bilayer vesicles consisting of, for example, 4000–16000 lipids Cooke et al. (2005), 4096 lipids Markvoort et al. (2009), and 5072 lipids Zheng et al. (2010). The curvature effect, arising from the nature of the bilayer, in these typical coarse-grained calculations is larger than that of GUVs with negligible curvature, and the vesicles prefer to flatten, creating a pore. Because of the curvature effect, the values for GUVs can vary compared with those reported here. However, the magnitude of the bending energy per lipid molecule of , estimated from typical bilayer bending stiffness and curvature radius , is much smaller than that of the attractive interaction potential of about and that of the electrostatic interaction potential of for . Therefore, the deviation from the values for GUVs should be small, and the obtained values here should be relevant to GUVs and biomembranes on the microscale.

The consistency of the topological changes observed in this study with our previous observations of characteristic pore formation in charged lipid vesicles Himeno et al. (2015) validates the plausibility of our coarse-grained MD simulation of charged lipid vesicles. The critical salt concentration for the charge-induced pore formation in charged lipid vesicles was determined to be on the order of , which is the typical salt concentration to which biomembranes are exposed. In other words, this salt concentration level is necessary to maintain the closed membrane morphology. Considering that biomembranes contain only negatively charged lipids and thus are subjected to electrostatic repulsion, the topological and morphological stabilities of biomembranes could also be physically maintained by the screened electrostatic interactions at a salt concentration of . Moreover, the morphological dynamics of biomembranes, such as exo- and endocytosis, may be finely controlled by the precise tuning of the electrostatic repulsive forces via the salt concentration near the critical point. This should be investigated by future studies, including the effects of lipid asymmetry between the two leaflets, other charged proteins, explicit solvents, and so on, with comparisons to experiments and more detailed simulations.

V Conclusion

In this study, we focused on the dynamics of binary charged lipid vesicles, which are typical soft materials physiologically relevant to the basic structure of biomembranes. The important physicochemical phenomena in this system are the lateral phase separation and the morphological changes. By systematically changing the interaction parameters and salt concentration, we found that the lateral phase separation is significantly delayed and/or inhibited in the presence of electrostatic repulsion. The analysis of domain size in neutral and charged lipid vesicles suggested a different mechanism of phase separation in neutral and charged vesicles, where the domain-growth exponents for spinodal decomposition and nucleation growth were observed, respectively. Moreover, the longer term calculation unveiled the morphological dynamics, such as pore formation and further transformation into disk, string, and bicelle shapes, depending on the attractive interaction among lipids and salt concentration, as summarized in the phase diagram. The MD visualization suggested that the pore formation in the charged-lipid-rich domain was initiated by local disturbances of the orientational order of bilayer membranes, which appears prior to the pore formation. Our findings provide fundamental insights into the structure, dynamics, and stability of lipid bilayer membranes and vesicles containing charged lipids, which are broadly seen in biomembranes, as well as artificial membranes mimicking the functional features of biomembranes.

Most of the calculations were performed using the parallel computer “SGI UV3000” provided by Research Center for Advanced Computing Infrastructure at JAIST. H.I. acknowledges support from Grants-in-Aid for the Japan Society for the Promotion of Science (JSPS) Research Fellowship for Young Scientists (JP13J01297) and Sasakawa Scientific Research Grant from The Japan Science Society (28-233). N.S. acknowledges support from the Grant-in-Aid for Scientific Reserch on Innovative Areas “Molecular Robotics” (JP15H00806) from the Ministry of Education, Culture, Sports, Science, and Technology of Japan (MEXT) and the Grant-in-Aid for Young Scientist (B) (JP26800222) from JSPS.


  • Eggeling et al. (2009) C. Eggeling, C. Ringemann, R. Medda, G. Schwarzmann, K. Sandhoff, S. Polyakova, V. N. Belov, B. Hein, C. von Middendorff, A. Schönle,  and S. W. Hell, Nature (London) 457, 1159 (2009).
  • Toulmay and Prinz (2013) A. Toulmay and W. A. Prinz, J. Cell Biol. 202, 35 (2013).
  • Simons and Ikonen (1997) K. Simons and E. Ikonen, Nature (London) 387, 569 (1997).
  • Simons and Sampaio (2011) K. Simons and J. L. Sampaio, Cold Spring Harbor Perspect. Biol. 3, a004697 (2011).
  • Suzuki et al. (2012) K. G. N. Suzuki, R. S. Kasai, K. M. Hirokawa, Y. L. Nemoto, M. Ishibashi, Y. Miwa, T. K. Fujiwara,  and A. Kusumi, Nat. Chem. Biol. 8, 774 (2012).
  • Veatch and Keller (2002) S. L. Veatch and S. L. Keller, Phys. Rev. Lett. 89, 268101 (2002).
  • Baumgart et al. (2003) T. Baumgart, S. T. Hess,  and W. W. Webb, Nature (London) 425, 821 (2003).
  • Komura and Andelman (2014) S. Komura and D. Andelman, Adv. Coll. Int. Sci. 208, 34 (2014).
  • Veatch and Keller (2003) S. L. Veatch and S. L. Keller, Biophys. J. 85, 3074 (2003).
  • Heberle and Feigenson (2011) F. A. Heberle and G. W. Feigenson, Cold Spring Harbor Perspect. Biol. 3, a004630 (2011).
  • Saeki et al. (2006) D. Saeki, T. Hamada,  and K. Yoshikawa, J. Phys. Soc. Jpn. 75, 013602 (2006).
  • Yanagisawa et al. (2007) M. Yanagisawa, M. Imai, T. Masui, S. Komura,  and T. Ohta, Biophys. J. 92, 115 (2007).
  • Stanich et al. (2013) C. A. Stanich, A. R. Honerkamp-Smith, G. G. Putzel, C. S. Warth, A. K. Lamprecht, P. Mandal, E. Mann, T. D. Hua,  and S. L. Keller, Biophys. J. 105, 444 (2013).
  • Vequi-Suplicy et al. (2010) C. C. Vequi-Suplicy, K. A. Riske, R. L. Knorr,  and R. Dimova, Biochim. Biophys. Acta 1798, 1338 (2010).
  • Shimokawa et al. (2010) N. Shimokawa, M. Hishida, H. Seto,  and K. Yoshikawa, Chem. Phys. Lett. 496, 59 (2010).
  • Blosser et al. (2013) M. C. Blosser, J. B. Starr, C. W. Turtle, J. Ashcraft,  and S. L. Keller, Biophys. J. 104, 2629 (2013).
  • Pataraia et al. (2014) S. Pataraia, Y. Liu, R. Lipowsky,  and R. Dimova, Biochim. Biophys. Acta 1838, 2036 (2014).
  • Himeno et al. (2014) H. Himeno, N. Shimokawa, S. Komura, D. Andelman, T. Hamada,  and M. Takagi, Soft Matter 10, 7959 (2014).
  • van Deenen et al. (1982) L. L. M. van Deenen, J. A. F. O. den Kamp, B. Roelofsen,  and K. W. A. Wirtz, Pure and Appl. Chem. 54, 2443 (1982).
  • Yeung et al. (2008) T. Yeung, G. E. Gilbert, J. Shi, J. Silvius, A. Kapus,  and S. Grinstein, Science 319, 210 (2008).
  • Dowhan (1997) W. Dowhan, Annu. Rev. Biochem. 66, 199 (1997).
  • Schlame (2008) M. Schlame, J. Lipid Res. 49, 1607 (2008).
  • May et al. (2002) S. May, D. Harries,  and A. Ben-Shaul, Phys. Rev. Lett. 89, 268102 (2002).
  • Shimokawa et al. (2011) N. Shimokawa, S. Komura,  and D. Andelman, Phys. Rev. E 84, 031919 (2011).
  • Okamoto et al. (2016) R. Okamoto, N. Shimokawa,  and S. Komura, Eur. Phys. Lett. 114, 28002 (2016).
  • Shimokawa et al. (2016) N. Shimokawa, H. Himeno, T. Hamada, M. Takagi, S. Komura,  and D. Andelman, J. Phys. Chem. B 120, 6358 (2016).
  • McMahon and Gallop (2005) H. T. McMahon and J. L. Gallop, Nature (London) 438, 590 (2005).
  • Beckers et al. (1989) C. J. Beckers, M. R. Block, B. S. Glick, J. E. Rothman,  and W. E. Balch, Nature (London) 339, 397 (1989).
  • Juhasz and Neufeld (2006) G. Juhasz and T. P. Neufeld, PLoS Biol. 4, e36 (2006).
  • Boroske et al. (1981) E. Boroske, M. Elwenspoek,  and W. Helfrich, Biophys. J. 34, 95 (1981).
  • Hamada et al. (2009) T. Hamada, Y. Hirabayashi, T. Ohta,  and M. Takagi, Phys. Rev. E 80, 051921 (2009).
  • Hamada et al. (2010) T. Hamada, R. Sugimoto, M. C. Vestergaard, T. Nagasaki,  and M. Takagi, J. Am. Chem. Soc. 132, 10528 (2010).
  • Khalifat et al. (2008) N. Khalifat, N. Puff, S. Bonneau, J.-B. Fournier,  and M. I. Angelova, Biophys. J. 95, 4924 (2008).
  • Saarikangas et al. (2009) J. Saarikangas, H. Zhao, A. Pykäläinen, P. Laurinmäki, P. K. Mattila, P. K. J. Kinnunen, S. J. Butcher,  and P. Lappalainen, Curr. Biol. 19, 95 (2009).
  • Ito et al. (2015) H. Ito, Y. Nishigami, S. Sonobe,  and M. Ichikawa, Phys. Rev. E 92, 062711 (2015).
  • Nishigami et al. (2016) Y. Nishigami, H. Ito, S. Sonobe,  and M. Ichikawa, Sci. Rep. 6, 18964 (2016).
  • Rajendran and Simons (2005) L. Rajendran and K. Simons, J. Cell Sci. 118, 1099 (2005).
  • Dall’Armi et al. (2013) C. Dall’Armi, K. A. Cevereaux,  and G. D. Paolo, Curr. Biol. 23, R33 (2013).
  • Sakuma et al. (2008) Y. Sakuma, M. Imai, M. Yanagisawa,  and S. Komura, Eur. Phys. J. E 25, 403 (2008).
  • Sakuma et al. (2010) Y. Sakuma, T. Taniguchi,  and M. Imai, Biophys. J. 99, 472 (2010).
  • Yanagisawa et al. (2008) M. Yanagisawa, M. Imai,  and T. Taniguchi, Phys. Rev. Lett. 100, 148102 (2008).
  • Hamada et al. (2007) T. Hamada, Y. Miura, K. Ishii, S. Araki, K. Yoshikawa, M. Vestergaard,  and M. Takagi, J. Phys. Chem. B 111, 10853 (2007).
  • Beales et al. (2011) P. A. Beales, C. L. Bergstrom, N. Geerts, J. T. Groves,  and T. K. Vanderlick, Langmuir 27, 6107 (2011).
  • Himeno et al. (2015) H. Himeno, H. Ito, Y. Higuchi, T. Hamada, N. Shimokawa,  and M. Takagi, Phys. Rev. E 92, 062713 (2015).
  • Markvoort et al. (2009) A. J. Markvoort, P. Spijker, A. F. Smeijers, K. Pieterse, R. A. van Santen,  and P. A. J. Hilbers, J. Phys. Chem. B 113, 8731 (2009).
  • Zheng et al. (2010) C. Zheng, P. Liu, J. Li,  and Y.-W. Zhang, Langmuir 26, 12659 (2010).
  • Venturoli et al. (2006) M. Venturoli, M. M. Sperotto, M. Kranenburg,  and B. Smit, Phys. Rep. 437, 1 (2006).
  • Noguchi (2009) H. Noguchi, J. Phys. Soc. Jpn. 78, 041007 (2009).
  • Cooke et al. (2005) I. R. Cooke, K. Kremer,  and M. Deserno, Phys. Rev. E 72, 011506 (2005).
  • Guigas and Weiss (2006) G. Guigas and M. Weiss, Biophys. J. 91, 2393 (2006).
  • Wang and Deserno (2016) X. Wang and M. Deserno, J. Phys. Chem. B 120, 6061 (2016).
  • Laradji and Kumar (2005) M. Laradji and P. B. S. Kumar, J. Chem. Phys. 123, 224902 (2005).
  • Lifshitz and Slyozov (1961) I. M. Lifshitz and V. V. Slyozov, J. Phys. Chem. Solids 19, 35 (1961).
  • Wagner (1961) C. Wagner, Z. Elektrochem. 65, 581 (1961).
  • Lipowsky and Dimova (2003) R. Lipowsky and R. Dimova, J. Phys.: Condens. Matter 15, S31 (2003).
  • Li et al. (2015) J. Li, H. Zhang, F. Qiu, Y. Yang,  and J. Z. Y. Chen, Soft Matter 11, 1788 (2015).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description