Solvent-free coarse-grained lipid model for large-scale simulations

Solvent-free coarse-grained lipid model for large-scale simulations


A coarse-grained molecular model, which consists of a spherical particle and an orientation vector, is proposed to simulate lipid membrane on a large length scale. The solvent is implicitly represented by an effective attractive interaction between particles. A bilayer structure is formed by orientation-dependent (tilt and bending) potentials. In this model, the membrane properties (bending rigidity, line tension of membrane edge, area compression modulus, lateral diffusion coefficient, and flip-flop rate) can be varied over broad ranges. The stability of the bilayer membrane is investigated via droplet-vesicle transition. The rupture of the bilayer and worm-like micelle formation can be induced by an increase in the spontaneous curvature of the monolayer membrane.

I Introduction

Amphiphilic molecules, such as lipids and detergents, self-assemble into various structures depending on the relative size of their hydrophilic parts: spherical or worm-like micelles, bilayer membranes, inverted hexagonal structures, and inverted micelles. Among these structures, the bilayer membrane of phospholipids has been intensively investigated, since it is the basic structure of the plasma membrane and the intracellular compartments of living cells, where the membranes are in a fluid phase and lipid molecules can diffuse in quasi-two-dimensional space. A vesicle (closed membrane) is considered to be a simple model of cells and it has applications in drug-delivery systems as a drug carrier.

Bilayer membranes exhibit many interesting phenomena such as shape deformation induced by phase separation or chemical reaction, membrane fusion, and membrane fission. The length scale of these phenomena varies from nm to m, since cells are m in diameter, whereas the thickness of a biomembrane is nm. To investigate the morphologies of cells and vesicles, the molecular structure is assumed to be negligible, and the bilayer membrane is described as a smoothly-curved mathematical surface Safran (1994); Lipowsky and Sackmann (1995); Seifert (1997). The information about the bilayer properties is only reflected in the values of the elastic parameters. To simulate the membrane with thermal fluctuations, a triangulated surface is widely used Gompper and Kroll (2004, 1997). An alternative model is a meshless membrane Drouffe et al. (1991); Noguchi and Gompper (2006a, b); Del Pópolo and Ballone (2008); Kohyama (2009); Liu et al. (2009); Füchslin et al. (2009); Yuan et al. (2010), where particles self-assemble into a membrane by anisotripic potential interactions. These models can reproduce m-scale dynamics of the bilayer membrane well but cannot treat a non-bilayer structure such as the stalk structure of a membrane fusion intermediate Jahn and Grubmüller (2002); Chernomordik and Kozlov (2008).

To simulate molecular-scale dynamics and the non-bilayer structure, a molecular model is required. Although computer technology has grown rapidly, the typical scale for recent simulations of the all-atom models is only ns dynamics of hundreds of lipid molecules. To simulate the membranes on longer and larger scales, various coarse-grained molecular models have been proposed (see review articles Müller et al. (2006); Venturoli et al. (2006); Klein and Shinoda (2008); Noguchi (2009); Marrink et al. (2009)). Recently, the potential parameters in some of the coarse-grained molecular models are tuned by atomistic simulations Marrink et al. (2004); Izvekov and Voth (2005); Arkhipov et al. (2008); Shinoda et al. (2008); Wang and Deserno (2010). In mapping of interaction parameters, one coarse-grained particle typically represents three or four heavy atoms and their accompanying hydrogen atoms. To further reduce the computational costs, larger segments (three or more segment particles per amphiphilic molecule) are employed, and the solvent is implicitly represented by an effective attractive potential between the hydrophobic segments Noguchi (2009); Noguchi and Takasu (2001a); Noguchi (2003); Brannigan et al. (2006); Deserno (2009); Farago and Grønbech-Jensen (2009). Model parameters are chosen to generate a bilayer membrane with reasonably realistic values of elastic properties. In this scale, it is difficult to take into account chemical details of lipids, such as an unsaturated bond in hydrocarbon chains. Instead, this type of models can be advantageous to capture the general features in the bilayer membrane, since the simplicity of the model can allow for wide ranges of variation of the membrane properties.

In this paper, we propose a solvent-free molecular model to pursue two purposes: (1) to represent the amphiphilic molecule in a size as small as possible and (2) to allow the variation in the membrane properties for wide ranges. A molecule consisting of many particles has a higher resolution than that with less particles but requires a smaller length unit and time step for simulations. Here, we consider a molecule that consists of a spherical particle and an orientation vector. It can reduce computational costs to simulate many molecules. In previous solvent-free models Noguchi (2009); Wang and Deserno (2010); Noguchi and Takasu (2001a); Noguchi (2003); Brannigan et al. (2006); Deserno (2009); Farago and Grønbech-Jensen (2009), the membrane properties are varied only in narrow ranges. On the other hand, in one of the meshless membrane models, the bending rigidity and the line tension of the membrane edge can be independently varied over wide ranges Noguchi and Gompper (2006a). This allows the conditions of vesicle formation and rupture to be controlled Noguchi and Gompper (2006b). In addition, it is easy to compare the simulation results with theoretical predictions. Such tuning capability is desired for molecular models.

In Sec. II, the lipid model and the simulation method are described. In Sec. III, the results and discussion are provided. The formation of a membrane and its stability for droplet-vesicle transitions are described in Sec. III.1. In Sec. III.2 and III.3, the calculation methods and the dependence of the static and dynamic properties on model parameters are described, respectively. The summary is given in Sec. IV.

Ii Model and Method

ii.1 Molecular model

Figure 1: (Color online) The cutoff function , the compact Gaussian weight function , and the repulsive potential .

In solvent-free lipid models, an amphiphilic molecule is typically represented by three or more particles Noguchi (2009). Here, we consider a lipid molecule with minimum () degrees of freedom for solvent-free molecular simulations. Each (-th) molecule has a spherical particle with an orientation vector , which represents a direction from the hydrophobic to the hydrophilic part. There are two points of interaction in the molecule: the center of a sphere and a hydrophilic point . The molecules interact with each other via the potential,


where , , , and is the thermal energy. The molecules have an excluded volume with a diameter via the repulsive potential,


with a cutoff at .

The second term in Eq. (1) represents the attractive interaction between the molecules. A multibody attractive potential is employed to mimic the “hydrophobic” interaction. This potential allows the formation of the fluid membrane over wide parameter ranges and fast lateral diffusion. Similar potentials have been applied in the previous membrane models Noguchi and Takasu (2001a); Noguchi and Gompper (2006a); Farago and Grønbech-Jensen (2009) and a coarse-grained protein model Takada et al. (1999). The potential is given by


where is chosen such that . The local particle density is approximately the number of particles in the sphere whose radius is .


where is a cutoff function,


with , , , and the cutoff radius (see Fig. 1). The potential acts as a pairwise attractive potential (, so that ) for and approaches a constant value () for . It is assumed that the hydrophobic parts have no contact with the implicit solvent (void space) at .

The third and fourth terms in Eq. (1) are discretized versions of tilt and bending potentials of the tilt model Hamm and Kozlov (1998, 2000), respectively. A smoothly truncated Gaussian function Noguchi and Gompper (2006a) is employed as a weight function


with , , and (see Fig. 1). All orders of derivatives of and are continuous at the cutoff radii. The weight is a function of (not ) to avoid the interaction between the molecules in the opposite monolayers of the bilayer. The average distance between the neighboring molecules in the same monolayer is , and the distance to the neighboring molecule in the other monolayer is (see Fig. 2(a)). Thus, these two potentials act between the neighboring molecules in the same monolayer but not between the monolayers. The tilt potential has the energy minimum in a completely flat membrane with no tilt deformation. Similar tilt potentials have been used in the meshless membrane models Drouffe et al. (1991); Del Pópolo and Ballone (2008); Kohyama (2009); Liu et al. (2009); Füchslin et al. (2009); Yuan et al. (2010). The same type of bending potential [the fourth term in Eq. (1)] was previously used to control the bending rigidity and the spontaneous curvature of the monolayer in the molecular simulations Noguchi (2003); Farago and Grønbech-Jensen (2009). Positive spontaneous curvature indicates that the hydrophilic head is larger than the hydrophobic tail of amphiphilic molecules. The bending rigidity is numerically calculated and compared with the estimation from the continuous description of the membrane in Sec. III.3.

Figure 2: (Color online) Probability distribution of (a) the positions and (b) the orientation of the molecules in the planar membrane at , , , , and . (a) The components of and are shown for . The dashed line represents of the molecules at . (b) The component of the molecular orientation is shown for , , , and . The error bars are displayed at several data points.

ii.2 Brownian dynamics

We simulated the membrane in the ensemble (constant number of molecules , volume , and temperature ) with periodic boundary conditions in a box with side length , , and . We employed Brownian dynamics (molecular dynamics with Langevin thermostat). The motions of the center of the mass and the orientation are given by underdamped Langevin equations:


where and are the mass and the moment of inertia of the molecule, respectively. The forces are given by and with the perpendicular component and a Lagrange multiplier to keep . According to the fluctuation-dissipation theorem, the friction coefficients , and the Gaussian white noises , obey the following relations: the average and the variance , where and . The Langevin equations are integrated by the leapfrog algorithm Allen and Tildesley (1987) with . First, the velocities are updated by




Then, the positions are updated by


We employed , , , , , , and the total number of the molecules to . The results are displayed with the length unit , the energy unit , and the time unit . The diffusion coefficient is normalized using the diffusion coefficient of an isolated molecule. The error bars of the data are estimated from the standard deviations of three to six independent runs.

Figure 3: (Color online) Sequential snapshots of the molecular self-assembly at , , , , , , and . (a) . (b) . (c) . (d) .
Figure 4: (Color online) Droplet-vesicle transition at , , , , and . The lower and upper lines represent the radius of gyration in increasing or decreasing, respectively. Sliced snapshots are also shown at .
Figure 5: (Color online) Formation of vesicles from a droplet at , , , , and . The tilt coefficient is gradually increased as . Sliced snapshots are shown at (a) (), (b) , (c) , (d) , and (e) (). All molecules are also shown for in (c’).

Iii Results and Discussion

iii.1 Self-assembly and membrane stability

Molecules self-assemble into spherical droplets at , i.e., when only the first two terms in Eq. (1) are taken into account. When the third term (the tilt potential) is added, the molecules can spontaneously form vesicles. Figure 3 shows the self-assembly of molecules from a random gas state. First, small clusters are formed; these clusters merge into disk-like micelles. Then, a large disk closes into a vesicle through a bowl-like shape (see Fig. 3(c)). Similar self-assembly processes have been observed in the previous simulations of molecular Noguchi and Takasu (2001a) and meshless Noguchi and Gompper (2006b) models.

In order to clarify the stability of three-dimensional aggregates and bilayer membranes, the morphologies of the aggregates are investigated as gradually increases or decreases. As increases, a spherical liquid droplet transforms into a bilayer vesicle. At the transition point (), the radius of gyration exhibits an abrupt increase as shown in Fig. 4. As decreases, the transition from a vesicle to a droplet occurs, however, the transition point () is much lower. Thus, a typical hysteresis for the first-order transition is observed. The rate of increase or decrease is sufficiently low (). We checked that the deviation of the transition points by the annealing rates is very small; between and with , , , , and . The transition points are not sensitive to the path of , since the difference of the results for and constant is smaller than their statistical errors.

Figure 6: (Color online) Shape transition points of molecular aggregates between the droplet and the bilayer membrane (vesicles and disks) for various (a) , (b) , (c) , and (d) . If not specified, , , , and . The solid and dashed lines represent data with increasing and decreasing , respectively. The open and filled symbols represent data with fixed and , respectively.
Figure 7: (Color online) Surface tension of a flat membrane at , , , , and . Dependence of on (a) the projected area per molecule and (b) the intrinsic area per molecule . The squares, triangles, and circles represent for , and , respectively. The error bars are smaller than the line thickness.

For large aggregates with , two vesicles or a vesicle with disks are formed instead of a single vesicle. Figure 5 shows an example of the formation of two vesicles. A void space is opened in the droplet, and a bilayer skirt is formed. Then, it is separated into two parts and forms two vesicles. If the separated membrane is small, a disk is formed. Since the larger droplets can have a clearer molecular layer on the surface, which prevents the shape change, higher is needed to trigger the shape transition [see Fig. 6(a)]. At , the coexistence region of the droplets and the vesicles is narrow, since the number of the molecules is not sufficient to form the surface and inside layers. The points of the droplet-bilayer transition are also dependent on , while they are almost independent of and . Thus, the bilayer stability is determined by the tilt and bending potentials but not by the attractive potential.

Let us discuss the condition required to form a stable bilayer. When all molecules have the same orientation , the bending potential energy becomes zero at . Thus, the bending potential with can have the minimum energy for any structure of the aggregate. Therefore, the spherical liquid droplet, which has the minimum surface area, would be the equilibrium state at instead of the bilayer. However, large vesicles with and planar membranes can maintain their bilayer structure as a metastable state even at with finite . Note that the capability to keep a pre-formed bilayer membrane does not guarantee the self-assembly to the bilayer. In particular, the periodic boundary condition is a strong constraint, which can keep the bilayer membrane as a thin liquid layer even at . In order to obtain the spontaneous formation of the bilayer membrane, is required.

iii.2 Calculation of membrane properties

To investigate the membrane properties, we formed a nearly planar membrane without edges or pores. The membrane area and the surface tension are varied by increasing or decreasing the projected area , where . The membrane has a clear bilayer structure (see Fig. 2). The intrinsic area of the tensionless membrane per molecule , the area compression modulus , and the half lifetime of the flip-flop motion are calculated from the flat membranes with . The bending rigidity and the diffusion coefficient are calculated at larger tensionless membranes with . The line tension of membrane edge is calculated from the strip of the flat membrane with .

Figure 8: (Color online) Spectra of undulation modes of nearly planar, tensionless membranes () at , , , , and . Results for calculated from the molecular positions () and from the averaged positions on a square mesh () are shown. The inset shows the dependence of on , which is used to extract the bending rigidity .
Figure 9: (Color online) Line tension of membrane edge at , , , , and . The circles represent calculated from a pore on the flat membrane at . The solid line represents calculated from the striped membrane at : .
Figure 10: (Color online) Time development of the probability difference of the molecules in the upper and lower monolayers at , , , , and . At the initial states (), and . The error bars are displayed at several data points.

The surface tension is given by Rowlinson and Widom (1982); Allen and Tildesley (1987)


with the diagonal components of the pressure tensor


where . When the potential interaction crosses the periodic boundary, the periodic image nearest to the other interacting molecules is employed. The intrinsic area of the membrane is larger than the projected area in the -plane due to the membrane undulations. We calculate from a square mesh with . The height of a mesh point is obtained from the weighted average of molecular position in the four neighbor cells, with and . Figure 7 shows the area dependence on the surface tension . The tension exhibits a roughly linear increase with the molecular area at . The compressed membrane with buckles out of plane and has the larger intrinsic area than the projected area . Similar dependence and buckling are obtained in the simulations of other molecular models den Otter (2005); Deserno (2009) and meshless models Noguchi and Gompper (2006a); Kohyama (2009).

The area of the tensionless membrane () is obtained by the minimization of , where the projected area is updated as every interval, where is the time average for . We use to and or . The area compression modulus is defined as


We calculate from the slope of - lines shown in Fig. 7(b).

The bending rigidity is calculated from the spectra of undulation modes of the planar membranes in Fourier space Safran (1994); Goetz et al. (1999); Lindahl and Edholm (2000),


Figure 8 clearly shows the dependence of the tensionless membrane. We calculate from the raw data (the particle position ), as well as from the square mesh with the same mesh-points which were used for the estimation of the intrinsic area . Averaging over the mesh removes most of the effects of the molecular protrusions. The bending rigidity is estimated from a fit of for , where the difference of two spectra is very small (see the inset of Fig. 8).

The line tension of the membrane edge is calculated from the strip of the flat membrane as Tolpekina et al. (2004); Reynwar and Deserno (2008)


since is the energy per unit length of the membrane edge, and the length of the membrane edge is . Since the striped membrane is tensionless, for solvent-free simulations. The tension and its error bar are estimated from the average and standard deviations for , , , and . Alternatively, can be also calculated from a circular pore on the flat membrane Tolpekina et al. (2004); Noguchi and Gompper (2006a). In this case, is balanced with the surface tension as . Since one has to estimate the pore radius , this method gives larger statistical errors as shown in Fig. 9. Therefore, we used the membrane strip for the calculation of .

Figure 11: (Color online) Parameter dependence of (a) the intrinsic area per molecule, (b) area compression modulus , (c) bending rigidity , (d) line tension , and (e) diffusion coefficient for the tensionless membrane at , , and . The circles and squares represent data for and , respectively.
Figure 12: (Color online) Parameter dependence of (a) , (b) , (c) , (d) , and (e) for , , and . The triangles, circles, and squares represent data for , , and , respectively.

In bilayer membranes, molecules can move laterally on a monolayer and transversely between upper and lower monolayers. The lateral diffusion coefficient of the molecules is calculated from the diffusion of the molecular projections in the plane; . In the fluid phase, the molecules exhibit a fast diffusion rate to .

The relaxation time of the transverse motion (flip-flop) between the upper and lower monolayers is measured from the relaxation of the labeled molecules Kornberg and McConnell (1971); Noguchi and Takasu (2001a). The differential equation of the probability () of the molecules, which belong to the upper (lower) monolayer, is given by , where . For planar membranes, , and at . When the molecules in the upper monolayer are initially labeled (), the probability decays as


with the half lifetime . Figure 10 shows that indeed follows the exponential decay in our simulations. Either the component of the position or orientation can be used to detect the nomolayer, to which a molecule belongs. Both and give the same probability distribution (entirely same for most of the parameters), since both the quantities have a clear minimum between two peaks (see Fig. 2).

Figure 13: (Color online) Bending rigidity dependence on (a) , (b) , and (c) for , , and . The solid lines with squares, circles, and triangles represent data for , , and , respectively. The dashed lines with crosses and diamonds represent data for and , respectively.
Figure 14: Parameter dependence of (a) , (b) , (c) , and (d) for , , and . The squares, circles, and triangles represent data for , , and , respectively.

iii.3 Parameter dependence of membrane properties

Figures 1118 show the parameter dependence of the properties of the tensionless membrane. The bilayer membrane is formed in the fluid phase over broad ranges of the parameters due to the multibody attractive potential. A gel phase is obtained only at a large value of in Eq. (3). The fluid-gel transition occurs at and for and , respectively [see jumps of in Fig. 11(e)]. As decreases, the lateral diffusion and flip-flop motion become faster, and the membrane elasticities (, , ) decrease [see Figs. 11 and 18(d)]. The intrinsic area per molecule ( molecules in one of monolayers) decreases with increasing and approaches the closest-packing area in the gel phase. In this paper, we focus on the fluid membrane and set , hereafter.

Figure 15: (Color online) Parameter dependence of (a) , (b) , (c) , and (d) for , , and . The squares and circles represent data for and , respectively.

The dependence on the strength of attraction in Eq. (1) is shown in Fig. 12. It has a tendency similar to dependence. The line tension can be varied by . At , is close to and the bilayer membrane is accompanied by free molecules (gas) with the average density of the gas . The molecules depart from the bilayer membrane and return. At with and , the bilayer membrane breaks and small micelles are formed ( molecules per micelle for ). On the other hand, no free molecules are seen at . The critical micelle concentration (CMC) of lipids is very low, and their chemical potential difference in solution and in membrane is typically more than per lipid Tanford (1980). Thus, the number of lipid molecules on vesicles is conserved in typical experiments. In order to keep the number of molecules on membrane constant during simulations, a sufficiently low CMC is required. We mainly use and , where the fluid membranes without free molecules are obtained.

Figure 16: (Color online) Line tension dependence on for (a) and (b) at and . (a) The diamonds, squares, triangles, and circles represent data for , and , respectively. (b) The squares, triangles, and circles represent data for , and , respectively.
Figure 17: (Color online) Parameter dependence of (a) , (b) , (c) , and (d) for , , and . The squares and circles represent data for and , respectively.

The bending rigidity of the bilayer membrane can be controlled by or . in Eq. (1). Figures 13(a) and (b) show the linear dependence of on and , respectively. When they are plotted together for , all lines are roughly overlapped on a line [see Fig. 13(c)]. A large deviation from the line is seen only at one data point at and (the leftmost triangle), where the bilayer structure is metastable. The bending rigidity is also weakly dependent on . The - curve maintains its shape and shifts upward with increasing as shown in Fig. 12(c). Thus, it is expressed by , where is an increasing function as , , and . The area compression modulus increases with increasing , while shows only slight dependence on for large (see Figs. 14 and 15). The other membrane properties , , and show weak dependence on and . Thus, can be varied by or without a large variation in . The modulus can be varied by .

Figure 18: (Color online) Half lifetime of flip-flop motion. (a) Dependence on at , , and . The dashed lines represent the free-energy barrier estimated by the orientation distribution shown in Fig. 2(b). (b) Dependence on at and . (c) Dependence on at , , and . (d) Dependence on at and .
Figure 19: (Color online) Sequential snapshots of vesicle rupture at , , , , , and . (a) . (b) . (c) .

The bending elastisity generated by the orientation-dependent potentials can be derived from the Helfrich theory for monolayer membranes. When the orientation vectors are equal to the normal vectors of the monolayers without tilt deformation, the bending and tilt energies are written by


in the continuum limit, where and are two principal curvatures of the monolayer. The first and second terms in Eq. (20) are the contributions of the bending and tilt potentials, respectively. The spontaneous curvature of the bending potential is given by . Noguchi (2003) The nearest-neighbor distance is obtained from the radial distribution function. By assuming the hexagonal packing of the molecules in the monolayers, the monolayer bending rigidities generated by the bending and tilt potentials are estimated as and , respectively. The bending rigidity of the monolayer is given by the sum of these . For the monolayer membrane, Eq. (20) gives the saddle-splay modulus and the spontaneous curvature with . Thus, the bending rigidity of the bilayer is estimated as from . This estimation of supports the linear dependence of the obtained simulation results on and . The prefactor () of gives the quantitative agreement, whereas the prefactor of is half of the numerical estimation. In the simulation, the thermal fluctuations induce molecular protrusion and tilt. These tilt fluctuations likely change the prefactor of to twice its value (). The attractive potential also adds a small bending resistance, .

The flip-flop time shows exponential dependence on the parameters of the tilt and bending potentials, while it has weak dependence on (see Fig. 18). The free-energy barrier between two monolayers would be the main factor to determine . It can be roughly estimated from the probability distribution of the molecular orientation as [see Fig. 2(b)]. The dashed lines in Fig. 18(a) represent calculated by this method. It gives very good agreements with , and thus, the flip-flop time is written by with . The barrier