Lipids with bulky head groups generate large membrane curvatures by small compositional asymmetries

Glycolipids such as GM1 have bulky head groups consisting of several monosaccharides. When these lipids are added to phospholipid bilayers, they generate large membrane curvatures even for small compositional asymmetries between the two leaflets of the bilayers. On the micrometer scale, these bilayer asymmetries lead to the spontaneous tubulation of giant vesicles as recently observed by optical microscopy. Here, we study these mixed membranes on the nanometer scale using coarse-grained molecular simulations. The membrane composition is defined by the mole fractions ϕ1 and ϕ2 of the large-head lipid in the two leaflets of the bilayer. Symmetric membranes are obtained for ϕ1 = ϕ2 ≡ ϕle, and asymmetric ones for ϕ1 ≠ ϕ2. In both cases, we compute the density and stress profiles across the membranes. The stress profiles are used to identify the tensionless states of the membranes. Symmetric and tensionless bilayers are found to be stable within the whole composition range 0 ≤ ϕle ≤ 1. For these symmetric bilayers, both the area compressibility modulus and the bending rigidity are found to vary non-monotonically with the leaflet mole fraction ϕle. For asymmetric bilayers, we compute the product of bending rigidity and spontaneous curvature from the first moment of the stress profile and determine the bending rigidities of the asymmetric membranes using the ϕle-dependent rigidities of the single leaflets. When we combine these results, the compositional asymmetry ϕ1 - ϕ2 is found to generate the spontaneous curvature (ϕ1 - ϕ2)/(0.63 ℓme) with the membrane thickness ℓme ≃ 4 nm. Therefore, the spontaneous curvature increases linearly with the compositional asymmetry. Furthermore, the small compositional asymmetry ϕ1 - ϕ2 = 0.04 leads to the large spontaneous curvature 1/(63 nm) and the increased asymmetry ϕ1 - ϕ2 = 0.2 generates the huge spontaneous curvature 1/(13 nm). These large values of the spontaneous curvature will facilitate future simulation studies of various membrane processes such as bud formation and nanoparticle engulfment.


I. INTRODUCTION
In spite of their diversity and molecular complexity, biomembranes have a universal architecture which is based on fluid bilayers of lipids and membrane proteins. 1 Within these bilayers, the hydrophilic head groups of the lipids are positioned between their hydrophobic chains and the aqueous solutions. The size of the lipid head groups can vary significantly: Phospholipids such as 1-palmitoyl-2-oleoyl-snglycero-3-phosphocholine (POPC) have a relatively small head (SH) group with a volume of about 0.5 cubic nanometers, while glycolipids such as the ganglioside GM1 have much bulkier head groups. 2 The glycolipid GM1, which has a head group consisting of four monosaccharides, is abundant in all mammalian neurons, 3 plays an important role in many neuronal processes and diseases, 4 and acts as a membrane anchor for various toxins, bacteria, and viruses. 5 Giant vesicles of POPC membranes doped with a small amount of GM1 have been recently studied experimentally a) URL: http://www.mpg.mpikg.de/th. b) lipowsky@mpikg.mpg.de and were observed to undergo spontaneous tubulation. 6,7 Such a spontaneous formation of membrane nanotubes provides direct evidence for bilayer asymmetry and spontaneous curvature. 8 This curvature can be generated by a variety of molecular mechanisms, including lipid-anchored macromolecules, 8 adsorption, 9 and depletion 10 of solutes and small molecules, membrane-bound proteins such as α-synuclein, 11 and Bin/Amphiphysin/Rvs-domain (BAR-domain) mimetics based on DNA origami. 12 Furthermore, all cell membranes have a compositional asymmetry in the sense that the two membrane leaflets differ in their composition. 13,14 Here, we consider such a compositional asymmetry for binary mixtures of lipid molecules with small and large head groups. The small-head (SH) lipids represent a phospholipid such as POPC, and the large-head (LH) lipids a glycolipid such as GM1. We model these lipids in a coarse-grained manner as displayed in Fig. 1, focusing on the different head group sizes and ignoring atomistic details. Therefore, we cannot expect to find quantitative agreement with experimental data for specific lipids. Nevertheless, our simulations show the same general trends as the experimental observations. Furthermore, when assembled into lipid bilayers, the model lipids in Fig. 1 FIG. 1. Two lipid species as used in the simulations: (a) Lipid SH has a small head group, consisting of three head group beads (H, red), and two hydrophobic chains, each of which is built up from six chain beads (C, green) and (b) Lipid LH has a large head group consisting of six H beads (blue) and two chains with six C beads. All H beads of both lipid species experience the same interactions with the other beads and likewise for all C beads. The different colors red and blue of the H beads are used to distinguish the two lipid species in simulation snapshots. generate large spontaneous curvatures for small compositional asymmetries and are thus useful to elucidate the influence of this curvature on various membrane processes such as bud formation and nanoparticle engulfment. We study these bilayer membranes by molecular simulations based on Dissipative Particle Dynamics (DPD). 15,16 This simulation method has been previously used to elucidate various aspects of bilayer membranes such as their interactions with nanoparticles, [17][18][19][20][21] membrane fusion, 22 receptor-mediated interactions between bilayers, 23 membrane poration by antimicrobial peptides, 24 and membrane curvature generated by Shiga toxin proteins. 25 In our DPD study, the membrane composition of the binary mixtures is described by the mole fractions φ 1 and φ 2 of the large-head lipids in the upper and lower leaflets of the bilayer. Symmetric membranes are obtained for mole fractions φ 1 = φ 2 ≡ φ le , and asymmetric ones for φ 1 φ 2 . For both symmetric and asymmetric membranes, we focus on tensionless membranes which we identify via their stress profiles. 26 For symmetric bilayers, we compute the average molecular area, the area compressibility modulus, and the bending rigidity as a function of mole fraction φ le within each leaflet. For asymmetric bilayers, we compute the product of bending rigidity and spontaneous curvature from the first moment of the stress profile and determine their bending rigidities from the φ le -dependent rigidities of the symmetric bilayers. Combining both results, we show that a relatively small bilayer asymmetry as described by a small value of |φ 1 − φ 2 | leads to a relatively large spontaneous curvature.
The paper is organized as follows. We first describe our molecular modeling and computational methods in Sec. II. Section II also defines the mole fractions φ 1 and φ 2 of the large-head lipids in the two leaflets of the bilayer membranes. In Sec. III, we discuss the density and stress profiles of symmetric and asymmetric bilayer membranes. In Sec. IV, we address the molecular areas of the lipid molecules. The elastic properties of symmetric bilayers are described in Sec. V and the spontaneous curvature of asymmetric bilayers in Sec. VI.

A. Molecular architecture and interactions
In lipid membranes, various quantities such as area compressibility, bilayer thickness, and bending rigidity are related to nanoscopic length and time scales. Here, we study these quantities using the coarse-grained DPD method 15 to simulate bilayer membranes. 16 In DPD, the particles are modeled as soft beads and their coordinates evolve according to Newton's equations of motion. All forces between the beads are pair-wise additive and conserve momentum which ensures a reliable description of hydrodynamics.
In this method, clusters of atoms are coarse-grained into beads with diameter, d, which provides the basic length scale as well as the force cutoff. We use three kinds of beads, water beads (W), hydrophilic head group beads (H), and hydrophobic chain beads (C) to construct the architecture of the model lipids as displayed in Fig. 1. We consider two types of lipids, SH and LH, that differ in the size of their head groups: the small-head (SH) lipids have a head group consisting of three H beads, while the head group of the large-head (LH) lipids is built up from six H beads. All H beads of both lipid species experience the same interactions with other beads. For simplicity, all beads were taken to have the same mass m 0 . The basic energy scale is provided by the thermal energy k B T, with Boltzmann's constant k B and temperature T, which implies the basic time scale τ = d 2 m 0 /(k B T ). The integration time step ∆t was taken to be ∆t = 0.01 τ as in previous studies. 9 Furthermore, the observed diffusion constant for lateral diffusion within the bilayer implies that the basic time scale τ is of the order of 1 ns.
In the DPD simulations, bead j exerts three pair-wise additive forces on bead i, a conservative force, a random force, and a dissipative force. The conservative force F co ij arises from both bonded and non-bonded interactions of the beads. The bonded interactions between adjacent beads in a lipid molecule are described by a Hookean spring, given by where r is the distance between the center-of-masses of two adjacent beads with the spring constant k r = 128 k B T /d 2 and equilibrium separation r eq = d/2. In addition, the two hydrophobic chains of the lipid molecule are stiffened by the bending potential where k φ = 15 k B T is the bending constant, θ is the tilt angle between two neighbouring bonds with equilibrium value θ eq = π corresponding to collinear bonds. In addition, all pairs of DPD beads experience conservative forces which are repulsive and described by where the force parameter f ij determines the effective hydrophobicity or hydrophilicity of the beads, r ij = r i − r j is the distance between the center-of-masses of the beads i and j, andr ij = (r i − r j )/r ij is the unit vector pointing from bead j to bead i. 27 The random and dissipative forces, F ra ij and F di ij , which provide the thermostat of the system, have the form and where γ ij is the friction coefficient, the variable ξ ij represents Gaussian white noise with ξ ij (t) = 0 and ξ ij (t)ξ i j (t ) = (δ ii δ jj + δ ij δ ji )δ(t − t ), and v ij = v i −v j is the relative velocity of particle i with respect to particle j. After combining all forces, we integrate Newton's equation of motion for each bead using a modified version of the velocity-Verlet algorithm 15 and study the time evolution of the system.

B. Leaflet compositions and mole fractions
In our simulations, we start from a pre-assembled bilayer containing a binary mixture of lipids with small and large head groups, SH and LH, see Fig. 1. Each small-head lipid has a head group consisting of three H beads whereas each largehead lipid has a head group with six H beads. Both lipid species have two chains with six C beads each. The bead-architecture of the small-head lipid as shown in Fig. 1(a) is identical with the bead-architecture of the lipids simulated in Ref. 9.
To focus on the bilayer asymmetry arising from the different compositions of the two leaflets of the membrane, both leaflets are taken to contain the same number of lipids but different mole fractions of the two lipids. As a consequence, both leaflets are characterized by the same average area per lipid as will be discussed in more detail in Sec. IV below. In this way, we separate the spontaneous curvature generated by the compositional asymmetry from the spontaneous curvature generated by different lipid densities in the two leaflets as studied previously in Ref. 9 for one-component lipid bilayers.
In most simulations, we used a total number of 29 2 = 841 lipid molecules in each leaflet. To distinguish the two leaflets, we introduce the Cartesian coordinate z perpendicular to the midplane of the bilayer. The location of this midplane is identified with the maximum of the C bead density corresponding to the hydrophobic core of the bilayer. The origin z = 0 of the z-coordinate is taken to be at this midplane. We then define the upper and the lower leaflet to be located at z > 0 and z < 0, respectively. Furthermore, quantities that refer to the upper leaflet will be indicated by the subscript 1, and those of the lower leaflet by the subscript 2.
The numbers of small-head and large-head lipids in the upper leaflet are denoted by N SH,1 and N LH,1 , and the corresponding numbers in the lower leaflet by N SH,2 and N LH,2 . As mentioned, we impose the constraint where N le represents the total number of lipids in each leaflet. For the simulations described below, this constraint on N le was conserved over the whole run time, which was typically of the order of 20 µs, because the two types of lipid molecules underwent essentially no flip-flops on these time scales. The compositional asymmetry is then defined by the two mole fractions and of the large-head group lipid LH in the upper and lower leaflets, respectively. In the following, we will study both symmetric and asymmetric bilayers. For symmetric bilayers, the two mole fractions φ 1 and φ 2 are equal and we will vary the leaflet mole fraction, For asymmetric bilayers, we will use the mole fraction φ 1 in the upper leaflet as the basic asymmetry parameter and vary this mole fraction for fixed φ 2 = 0. In the latter case, we will explore the whole range of possible asymmetries and simulate bilayers with mole fractions φ 1 = 0.024, 0.048, 0.072, 0.096, 0.2, 0.4, 0.6, 0.8, 1. The case φ 1 = 0 was previously studied in Ref. 9.

C. Simulation parameters
The interaction strength, f ij , and the friction coefficient, γ ij , assigned to different pairs of beads are given in Table I  and Table II, respectively. These parameter values are chosen as in previous studies 9,22 in order to obtain a well-ordered bilayer in its fluid state with a bending rigidity κ 15 k B T. As mentioned, all H beads of both lipid species experience the same interactions with other beads, i.e., the force parameters f HH , f HC = f CH , and f HW = f WH apply to all H beads of both SH and LH lipids.  Using the simulation parameters as described above, we performed DPD simulations in the NVT ensemble. The standard DPD parametrization of water corresponds to the bulk density ρ W = 3/d 3 of the water beads, which matches the water compressibility at room temperature T = 298K. 15 To suppress finite size effects as observed in Ref. 28, the bead density away from the bilayer is always chosen to be equal to the bulk water density ρ W = 3/d 3 .
The initial state of the simulations consists of a preassembled lipid bilayer parallel to the xy-plane with the z-axis being perpendicular to this plane. The initial midplane of the bilayer is located at z = 0. For a bilayer consisting of only small-head lipids, the lateral size L x = L y ≡ L and the height L z of the simulation box were chosen to be L = 32 d and L z = 48 d. 9 The inclusion of large-head lipids requires an increased area per lipid to obtain a tensionless membrane. Depending on the mole fraction φ 1 of LH lipids, the lateral size of the simulation box varied within the interval 32d ≤ L ≤ 35.8 d, while the height of the box was L z 48 d, with small variations arising from the two constraints that both leaflets contain the same number N le of lipids and that the bead density away from the bilayer has the fixed value 3/d 3 .

D. Stress profiles and tensionless membranes
To determine the mechanical tension within the membrane, we calculate the stress profile 22,26 from the diagonal elements Σ T (z) = Σ xx (z) = Σ yy (z) and Σ N = Σ zz of the stress tensor, which describe the local stress components tangential and normal to the midplane of the bilayer, respectively. The corresponding components P T (z) and P N of the pressure tensor have the same magnitude as Σ T (z) and Σ N but the opposite sign. The off-diagonal elements of the pressure tensor are zero since the bilayer membrane is in a fluid state and does not sustain any shear deformations. The mechanical tension Σ is then obtained by integrating the stress profile over z, i.e., by To obtain the stress profile s(z) by simulations, we apply the method of Goetz and Lipowsky. 26 According to this method, the bilayer membrane and the corresponding stress tensor can be considered to be isotropic in all the directions parallel to the bilayer midplane, i.e., to the xy-plane. Hence we divide the simulation box into thin slices parallel to the xy-plane, and the stress tensor is averaged over each of these slices using Heaviside step functions. The mechanical tension Σ can be decomposed according to with the leaflet tensions defined by 9 and These integral expressions depend on the location of the midplane at z = 0. For both symmetric and asymmetric bilayers, we define this location by the maximum of the C bead density which corresponds to the hydrophobic core of the bilayer. Note that both leaflet tensions contain contributions from bead-bead interactions between lipids that belong to different leaflets. In principle, one could also decompose the mechanical tension Σ into more than two contributions, arising (i) from the interactions between the H and C beads of only those lipids that belong to the upper leaflet, (ii) from the bead-bead interactions between the lipids in the lower leaflet, (iii) from the bead-bead interactions between lipids in different leaflets, and (iv) from the interactions between the lipids and the water beads. In the present study, we did not pursue such a more elaborate decomposition of the mechanical tension because we are primarily interested in the fluid-elastic properties of the bilayers and these properties do not depend on how we decompose the mechanical tension.
Using the leaflet tensions as defined in Eqs. (13) and (14), a symmetric bilayer is characterized by identical leaflet tensions which we denote by Furthermore, a tensionless membrane is defined by vanishing mechanical tension, i.e., by For symmetric and tensionless bilayers, this condition is equivalent to tensionless leaflets with Σ le = Σ 1 = Σ 2 = 0 (symmetric and tensionless).
For asymmetric membranes, tensionless leaflets imply a tensionless bilayer but a tensionless bilayer does not, in general, imply tensionless leaflets as will be demonstrated further below; see Fig. 5. All symmetric and asymmetric bilayer membranes discussed in the following are tensionless with Σ = 0 unless we explicitly mention a nonzero value of Σ.
In general, one might also consider bilayers for which both leaflet tensions Σ 1 and Σ 2 vanish separately. However, for the compositional asymmetry considered here, it is not possible to replace the constraint Σ 1 + Σ 2 = 0 in Eq. (18) by the two constraints Σ 1 = 0 and Σ 2 = 0. Indeed, in order to separate the bilayer asymmetry induced by different leaflet densities 9 from the bilayer asymmetry arising from different leaflet compositions, we impose the constraint that both leaflets contain the same number N le of lipid molecules, see Eq. (6), which implies that the overall lipid density has the same value in both leaflets. The latter constraint is, in general, incompatible with the two constraints Σ 1 = 0 and Σ 2 = 0; see Fig. 5 below.

III. DENSITY AND STRESS PROFILES
A. Symmetric bilayer membranes

Density profiles of symmetric bilayers
First, we discuss the density profiles of tensionless and symmetric bilayers as displayed in the upper panels of Fig. 2 for mole fractions φ le = 0, 0.4, and 1. The profiles are plotted as a function of the Cartesian coordinate z perpendicular to the bilayer midplane at z = 0. For visual clarity, we combine the densities of all chain beads into a single density profile labeled by T, but distinguish the head bead densities H SH and H LH corresponding to the two lipid species with small and large head groups. The density profiles show that the lipids form proper bilayer structures, with the water beads being completely excluded from the hydrophobic core of the bilayer. The head beads are located at the interface between the chain and water beads. The two head group layers are well separated, with essentially no overlap. Away from the bilayer, the water density profile ρ W (dashed line, turquoise) approaches the bulk density ρ = 3/d 3 which is kept fixed in all simulations. In the case of symmetric membranes with only small-head lipids, corresponding to φ le = 0, we obtain the same density profiles for the head and chain beads as in Ref. 9.
As we increase the mole fraction φ le of the LH lipids with the large head group, both the density profile of the water beads and the combined density profile of the chain beads hardly change. Likewise, the density profile of the H SH beads for φ le = 0 is rather similar to the density profile of the H LH beads for φ le = 1. At intermediate φ le -values, the two peaks of the H SH density profile are somewhat closer to the bilayer midplane compared to the two peaks of the H LH density profile. Thus, the density profiles undergo only small changes as we increase the mole fraction φ le . By contrast, the stress profiles change quite significantly.

Stress profiles of symmetric bilayers
In the lower panels of Fig. 2, we display the stress profiles as calculated for symmetric bilayers with equal mole fractions φ 1 = φ 2 = φ le in both leaflets. The three panels correspond to the values φ le = 0, 0.4, and 1.0. In all cases, the bulk water region is homogeneous and does not contribute to the stress profile. For φ le = 0 and 0.4, see lower panels of Figs. 2(a) and 2(b); the hydrophobic core regions are characterized by a negative stress or positive pressure between the chain beads. As we increase the number of LH lipids, the mean area per lipid increases and the pressure between the chain beads decreases. For φ le = 1, see Fig. 2(c), the bottom panel; the stress profile of the tensionless membrane is even inverted with a positive stress or negative pressure within the hydrophobic core and a negative stress or positive pressure within the head group layers. Note that the stresses within the hydrophobic core and the head group layers must cancel each other in order to satisfy the condition that the mechanical tension Σ = 0 as in Eq. (16).
For φ le = 0, corresponding to only small-head lipids, the stress profile exhibits outer peaks at the interfaces between the head and water beads as well as a deep minimum in the hydrophobic core of the bilayer; see lower panel of Fig. 2 The symmetry implies that both leaflets contain the same mole fraction φ le of the LH lipids with a large head group: (a) Only SH lipids with a small head group corresponding to φ le = 0, (b) A mixture of SH and LH lipids with φ le = 0.4, and (c) Only LH lipids corresponding to φ le = 1. The top panels display the density profiles. The W bead density (dashed line, turquoise) drops to zero in the hydrophobic core of the bilayer where the C bead density (dotted line, green) has a pronounced maximum which defines the midplane of the bilayer. The head beads H SH of the SH-lipids (continuous line, red) and the head beads H LH of the LH-lipids (dashed-dotted line, dark blue) are located at the interface between the chain and the water beads. The three bottom panels display the corresponding stress profiles. head groups in contact with the water are compressed. In the latter case, the stress profile has still two peaks with a positive stress value at the two interfaces between the chain and head beads. Thus, for φ le = 0, the head group layers are stretched and would like to shrink, while the hydrophobic core is compressed and would like to expand. By contrast, for φ le = 1, the head group layers are compressed and would like to expand, whereas the hydrophobic core is stretched and would like to shrink.
The stress profile for φ le = 1 as shown in Fig. 2(c) is somewhat unusual, and one may wonder about the stability of the bilayer structure. Indeed, the large-head lipids could prefer to assemble into cylindrical micelles as has been observed in simulations of amphiphilic molecules with a single and flexible chain. 26 However, the preassembled membranes with φ le = 1 as studied here kept their bilayer structure during the whole duration of our simulations which had a typical run time of 20 µs.

Density profiles of asymmetric bilayers
In addition to the density profiles of symmetric membranes, we determined the density profiles for compositionally asymmetric membranes, by varying the mole fraction φ 1 of the large-head lipids in the upper leaflet for fixed φ 2 = 0 in the lower leaflet. The corresponding density profiles are shown in the upper panels of Fig. 3. Inspection of these panels reveals that the total density of the beads in the upper leaflet with z > 0 increases with increasing mole fraction of the large-head lipids. This increase is a direct consequence of our procedure to change the mole fraction φ 1 in the upper leaflet by replacing small-head by large-head lipids, keeping the total lipid number N SH,1 + N LH,1 in the upper leaflet fixed and equal to the total lipid number N SH,2 + N LH,2 in the lower leaflet.

Stress profiles of asymmetric bilayers
Next, we describe the stress profiles of asymmetric bilayers which were studied for mole fractions φ 1 = 0.096, 0.4, and 0.8 in the upper leaflet keeping φ 2 = 0 in the lower leaflet. The corresponding stress profiles are shown in the lower panels of Fig. 3. Comparing the lower and upper panels of this figure, we find that an increase in φ 1 with fixed φ 2 = 0 leads to an increase of the combined density of head and chain beads within the upper leaflet and a concomitant decrease of the positive stress or tension within this leaflet. These changes can be interpreted in terms of an increased repulsion between the head beads in the upper leaflet until this repulsion dominates and the stress in the upper leaflet becomes negative, corresponding to a positive pressure between the head beads.
The larger the mole fraction of the large-head lipids in the upper leaflet, the smaller the associated leaflet tension Σ 1 ; see lower panels of Fig. 3. As a consequence, the membrane has the tendency to bend (or bulge) toward the upper leaflet which implies a positive preferred or spontaneous curvature as discussed in Sec. VI below.

A. Molecular areas within individual leaflets
For each leaflet, the molecular areas A SH and A LH of the two lipid species depend on the mole fraction and the mechanical tension within this leaflet. Thus, we have with k = 1 for the upper leaflet and k = 2 for the lower leaflet. The total leaflet area A k is then given by and the average molecular area A k by with k = 1, 2.
Because the bilayer membrane spans the simulation box, the projected areas A 1 and A 2 of the upper and lower leaflets are equal with A 1 = A 2 = L x L y = L 2 . Likewise, the average molecular areas A 1 and A 2 have identical values, but have different decompositions into contributions from the two lipid species SH and LH for different mole fractions φ 1 φ 2 and/or different leaflet tensions Σ 1 Σ 2 as follows from Eq. (21).

B. Symmetric and tensionless bilayers
For symmetric bilayers, both leaflets contain the same numbers N SH and N LH of SH and LH lipids, which implies identical mole fractions φ 1 = φ 2 = φ le , and both leaflets experience the same leaflet tension Σ 1 = Σ 2 = Σ le as in Eq. (15). Furthermore, symmetric and tensionless membranes are characterized by Σ 1 = Σ 2 = Σ le = 0 as in Eq. (17). For such membranes, the expression in Eq. (21) leads to the average molecular area in both leaflets. The functional dependence of this molecular area on the mole fraction φ le is displayed in Fig. 4. For φ le = 0 and φ le = 1, we can directly measure the limiting values A 0 (φ le = 0) = A SH (0, 0) 1.217 (23) and of the two lipid species. If the molecular areas of the two lipid species were independent of the composition, we would obtain the average molecular areâ which is linear in the mole fraction φ le ; see Fig. 4. Comparison of the average molecular areas A 0 andÂ 0 shows that A 0 is reduced, for intermediate values of φ le , compared toÂ 0 .

C. Asymmetric and tensionless bilayers
For asymmetric bilayers, the upper and lower leaflet differ in their composition as described by different mole fractions By definition, the midplane is located at the maximum of the C bead density. Because the upper leaflet is compressed and the lower leaflet is stretched, the membrane prefers to bend (or bulge) toward the upper leaflet. φ 1 and φ 2 . Likewise, the leaflet tensions Σ 1 and Σ 2 are also different. Furthermore, asymmetric and tensionless membranes are characterized by Σ 2 = −Σ 1 as in Eq. (18). Therefore, asymmetric and tensionless membranes are characterized by the leaflet tension Σ 1 = Σ 1 (φ 1 , φ 2 ). Furthermore, for the special case with φ 2 = 0 as considered here, the leaflet tension Σ 1 is determined by the mole fraction φ 1 alone as displayed in Fig. 5.
It then follows from the general expression in Eq. (21) that the average molecular area A 1 of a lipid in the upper leaflet is given by (26) and the average molecular area A 2 of a lipid in the lower leaflet by (27) Note that the molecular area A 1 depends on the leaflet tension Σ 1 and the molecular area A 2 on the leaflet tension Σ 2 = −Σ 1 . Thus, if the upper leaflet is compressed with Σ 1 < 0 as in Fig. 5, the lower leaflet is stretched with Σ 2 = −Σ 1 > 0 and vice versa.

V. ELASTIC PROPERTIES OF SYMMETRIC BILAYERS
Two important elastic properties of bilayer membranes are provided by their area compressibility K A and by their bending resistance or bending rigidity κ. Both quantities satisfy a simple linear relationship that also involves the membrane thickness me .

A. Area compressibility modulus
The area compressibility modulus K A of a symmetric bilayer determines the fractional change in the projected area per lipid, A sy , from its optimal value A sy = A 0 at constant temperature when the membrane experiences the mechanical tension Σ. Here and below, the subscript "sy" stands for "symmetric." Thus, the elastic modulus K A is defined via the linear relationship which implies To obtain the area compressibility modulus from simulation data, we considered symmetric membranes with different mole fractions φ le = φ 1 = φ 2 and determined the mechanical tension Σ as a function of the projected lipid area A sy . Two examples for the functional dependence of Σ on A sy are displayed in Fig. 6. The area compressibility modulus K A then follows from the slope of the function Σ = Σ(A sy ) at the tensionless state with A sy = A 0 . To determine this slope, we used the data points with molecular areas A sy that belonged to the range (A sy − A 0 )/A 0 < 0.1. The area compressibility modulus K A depends on the membrane composition as described by the mole fraction φ le of the large-head lipids; see Fig. 8(a) and fourth row of Table III. Inspection of Fig. 8(a) shows that K A varies with φ le in a nonmonotonic manner. For φ le = 0 corresponding to small-head lipids only, we obtain K A 27.1 ± 0.4 k B T /d 2 . For small values of φ le , the area compressibility decreases until it reaches a minimum at φ le 0.1 with K A 25.4 ± 0.2 k B T /d 2 . A further increase in φ le leads to an increase of K A up to K A 34.2 ± 0.1 k B T /d 2 for φ le = 1. Thus, as we increase the mole fraction φ le of the large-head lipids, the bilayer softens for small values of φ le , but stiffens for larger φ le -values.

B. Thickness of bilayer membranes
We define the thickness me of the bilayer membrane by the distance between the two peaks of the density profile ρ H for the head beads. In the presence of two different head groups, we take the combined density profile of both head groups. The values for the membrane thickness me obtained in this manner are displayed in Fig. 8(b) and in the fifth row of Table III. For a symmetric membrane with φ le = 0, corresponding to only small-head lipids, we find me 5 d. The membrane thickness remains essentially constant up to φ le 0.6 and then decreases until it reaches the value me 4.8 d for φ le = 1.
The decrease of the membrane thickness for larger mole fractions of the large-head lipids is consistent with the increase in the average molecular area A 0 of the lipids as described in Sec. IV B and displayed in Fig. 4. For real lipid bilayers with a single phospholipid component, the separation of the two head group layers is of the order of 4 nm which implies that the bead diameter d 0.8 nm in physical units.

C. Bending rigidity of symmetric bilayers
For symmetric bilayers, the bending rigidity κ = κ sy , which represents one key parameter of curvature elasticity, 29,30 was determined in two ways. First, we used the simple relationship Inserting the measured values of K A and me into this expression, we obtain the values for the bending rigidity κ sy as displayed in Fig. 8(c) and in the penultimate row of Table III, with an uncertainty of less than 3 percent. Furthermore, we find that the bending rigidity obtained in this manner varies in a non-monotonic manner as a function of the mole fraction φ le , reflecting the non-monotonic dependence of the area compressibility modulus on φ le . The numerical coefficient 1/48 in Eq. (30) was derived in Ref. 31 using classical elasticity theory for thin solid-like films and considering the limit in which the two-dimensional shear modulus vanishes. Using a polymer brush model, Rawicz et al. 32 obtained the same relationship as in Eq. (30), but with the numerical coefficient 1/48 replaced by 1/24. However, the simulation study of a one-component bilayer 31 confirmed the value 1/48 within an accuracy of about 10 percent by measuring the three parameters κ sy , K A , and me independently.
To corroborate the accuracy of the relationship in Eq. (30) for binary SH/LH bilayers with different values of the mole fraction φ le , we also studied the thermally excited shape fluctuations of the membranes, using the Fourier mode analysis introduced in Ref. 31. The fluctuations were decomposed into a discrete set of Fourier modes with the wavenumber q = (q x , q y ) = (2π/L ) (n x , n y ) and integers n x and n y . The corresponding fluctuation spectra S(q) were determined as a function of q ≡ |q| = q 2 x + q 2 y for symmetric and tensionless membranes with different mole fractions φ le ; see Fig. 7. Curvature elasticity predicts that these fluctuation spectra behave as 9,31 S(q) ≈ k B T κ sy q 4 for small q and Σ 0, and thus depend only on a single parameter, the dimensionless bending rigidity κ sy /(k B T ). Fitting the simulated Fourier mode spectra S(q) in Fig. 7 to Eq. (31), we obtain the κ sy -values corresponding to the solid lines (red) in Fig. 7 and displayed in the last row of Table III, with an uncertainty of less than 10 percent. In Fig. 7, the broken lines (black) represent the functional form of S(q) as given by Eq. (31) using the κ sy -values computed via Eq. (30); see the penultimate row of Table III. Inspection of Fig. 7 reveals that the broken line is almost on top of the solid line for each mole fraction φ le , directly demonstrating the good agreement between the results of the two computational methods. Furthermore, comparing the κ syvalues in the last two rows of Table III leads to the conclusion that these two values differ by less than 10 percent for all four mole fractions for which the Fourier mode analysis has been performed. As a consequence, the results of our Fourier mode analysis are again consistent with the prefactor 1/48 in Eq. (30), but exclude the prefactor 1/24 as proposed in Ref. 32. In the following, we will focus on the κ sy -values computed from Eq.  Table III. For each mole fraction φ le , the solid line is almost on top of the broken one, demonstrating the good agreement between the two computational methods used to determine the bending rigidity κ sy . smaller uncertainty than those deduced from the Fourier mode analysis.

D. Bending rigidities of bilayer leaflets
When the bilayer membrane is bent, both leaflets are bent simultaneously. We now denote the bending rigidities of the upper and lower leaflets by κ 1 and κ 2 , respectively. The composition-dependent bending rigidity κ sy (φ le ) of the symmetric bilayer is then given by the sum of the two leaflet rigidities, i.e., κ sy (φ le ) = κ 1 (φ le ) + κ 2 (φ le ) (32) with κ 1 (φ le ) = κ 2 (φ le ) = 1 2 κ sy (φ le ) (symmetric bilayer). (33) Thus, using the values for κ sy (φ le ) for symmetric bilayers, we can directly obtain the values κ 1 (φ le ) and κ 2 (φ le ) for the bending rigidities of the two leaflets. This decomposition will be used further below to estimate the bending rigidity for asymmetric bilayers as well.

E. Composition dependence of bending rigidity
As displayed in Fig. 8(c) and in the penultimate row of Table III, the bending rigidity has the value κ sy = 14.1 ± 0.2 k B T or 5.8 ± 0.08 × 10 −20 J at room temperature for φ le = 0, i.e., for a one-component membrane consisting of small-head lipids only. As we increase the mole fraction φ le of the large-head lipids to φ le = 0.024, the bending rigidity FIG. 8. Symmetric and tensionless membranes: (a) Area compressibility modulus K A , (b) Membrane thickness me , and (c) Bending rigidity κ sy versus mole fraction φ le of the large-head lipids in the two bilayer leaflets. In panel c, the full data points (red) were obtained via Eq. (30), and the open data points (blue) from a fit to Eq. (31). Both K A and κ sy vary with φ le in a nonmonotonic manner. The bilayer thickness me decreases for large values of φ le 0.6, reflecting the increased area per lipid, see Fig. 4, and the concomitant reduction of the end-to-end distance of the lipid chains. The broken lines represent linear interpolations between neighboring data points. This interpolation procedure produces the apparent kinks at φ le = 0.6, which should smoothen out when we include more data points in the vicinity of this φ le -value. becomes smaller and attains the value κ sy = 13.95 ± 0.3 k B T or 5.73 ± 0.12 × 10 −20 J; see Fig. 8(c). Increasing the mole fraction to φ le = 0.096, the bending rigidity is decreased to κ = 13.66 ± 0.11 k B T or 5.6 ± 0.04 × 10 −20 J. As we increase the mole fraction to values above φ le 0.2, the bending rigidity increases again and reaches its maximal value κ = 15.25 ± 0.1 k B T or 6.3 ± 0.04 × 10 −20 J for φ le = 1, corresponding to large-head lipids only.
Increasing the mole fraction from φ le = 0 to φ le = 0.048, the bending rigidity obtained from our simulations decreases by about 5 percent; see Fig. 8(c) and penultimate row of Table III. The same trend of the bending rigidity for small mole fractions has been observed in an experimental study of POPC/GM1 bilayers by fluctuation analysis and micropipette aspiration. 33 The uncertainty of the rigidity values determined by these two experimental methods was, however, quite large. More recently, POPC/GM1 bilayers were also studied by tube pulling experiments 7 which led to the conclusion that the bending rigidity remains essentially constant in the small mole fraction regime and has a value of about 20 k B T with an uncertainty of 10 percent. The latter value has the same order of magnitude as the κ sy -values deduced from our simulations for 0 ≤ φ le ≤ 0.048; see penultimate row of Table III.

A. Negative first moment of stress profile
An asymmetric bilayer membrane prefers to bend in a certain manner as described by the preferred or spontaneous curvature. For a tensionless membrane with Σ = 0 as in (16), the spontaneous curvature can be deduced from the first moment of the stress profile s(z) using the relation 9,10,34 Thus, the negative first moment of the stress profile is proportional to the product of bending rigidity κ and spontaneous curvature m. A symmetric bilayer has a symmetric stress profile with s(−z) = s(z) which implies that the integral on the right-hand side of (34) vanishes, whereas an asymmetric stress profile with s(−z) s(z) usually leads to a nonzero integral. The relation (34) has been previously used to compute the parameter combination 2κm for planar bilayers with asymmetric adsorption 9 and depletion 10 layers. In the latter case, an analytical solution for hard-core interactions showed that the spontaneous curvature can indeed be obtained from the properties of a planer bilayer. One should also note that the relationship in Eq. (34) does not involve any additional assumptions about the individual leaflets. In particular, this relationship does not depend on the leaflet tensions Σ 1 and Σ 2 as defined in Eqs. (13) and (14). Likewise, it does not involve any assumptions about the pivotal surfaces of the two leaflets which have been studied in Ref. 35.
For the binary mixtures considered here, we obtain the parameter combination 2κm as plotted in Fig. 9(a)  upper leaflet and no large-head lipids in the lower leaflet; see also the numerical values of 2κm in Table IV. Inspection of Fig. 9(a) shows that the parameter combination 2κm increases monotonically and almost linearly with the mole fraction φ 1 . In order to obtain the spontaneous curvature m from these data for 2κm, we now need to determine the bending rigidity κ for asymmetric membranes.

B. Bending rigidity of asymmetric membrane
For symmetric membranes, the bending rigidity can be decomposed into two contributions from the individual leaflets as in (32). We now generalize this decomposition to asymmetric membranes which leads to κ ≡ κ 1 (φ 1 ) + κ 2 (φ 2 ) . (35) This expression can be evaluated using κ 1 (x) = 1 2 κ sy (x) and κ 2 (x) = 1 2 κ sy (x) as in (33). We then obtain with κ sy = κ sy (φ le ) as obtained for symmetric bilayers; see Fig. 8(c) and the numerical values in Table III. Using the expression (36), we obtained the κ-values for asymmetric membranes with φ 1 > 0 and φ 2 = 0 as plotted in Fig. 9(b), with the numerical values given in Table IV.

C. Composition dependence of spontaneous curvature
The values of the spontaneous curvature m as obtained from the parameter combination 2κm and the bending rigidity κ of the asymmetric membranes are displayed in Table IV and plotted in Fig. 9(c) as open diamonds (blue). For comparison, Table IV also includes the values of the estimate m for the spontaneous curvature as obtained from the parameter combination 2κm with κ replaced by κ sy for symmetric bilayers. The m -values are plotted in Fig. 9(c) as filled circles (red). The relative difference |m − m|/m is always smaller than 10 percent which provides an estimate for the accuracy of the m-values obtained here.
The m-values in Fig. 9(c) can be well fitted by the linear relation where the second equality follows from the membrane thickness me = 5 d for φ 1 ≤ 0.6 as displayed in Fig. 8(b). The positive value of m implies that the bilayer prefers to bend (or bulge) toward the upper leaflet as one would expect intuitively because the large-head lipids in the upper leaflet want to occupy more space. The simulations have been performed for the special case φ 1 > 0 and φ 2 = 0. By swapping the indices 1 and 2, we directly TABLE IV. Asymmetric membranes with mole fraction φ 1 of the large-head lipids in the upper leaflet and φ 2 = 0 in the lower leaflet: Numerical values of the leaflet tension Σ 1 , the parameter combination 2κm, the bending rigidity κ as determined via Eq. (36), the spontaneous curvature m as obtained from the values for 2κm and κ, the estimate m = 2κm/κ sy for the spontaneous curvature with κ sy from Table III, and the relative difference δm ≡ |m − m |/m which is always smaller than 8 percent. The parameter values displayed here are also used for the plots in Fig. 9.
The spontaneous curvature as given by (39) can be quite large. For small compositional asymmetries with φ 1 − φ 2 = 0.01 and φ 1 − φ 2 = 0.04, for example, we obtain m = 1/(252 nm) and m = 1/(63 nm), respectively, where we used the membrane thickness me = 4 nm. The latter value is quite significant and would lead to cylindrical nanotubes with a diameter of only 63 nm. 36 Furthermore, for the mole fraction difference φ 1 − φ 2 = 0.4, Eq. (39) leads to m = 1/(1.58 me ) which is of the order of the inverse membrane thickness. For even larger values of φ 1 − φ 2 and thus larger values of m, the bilayer is likely to become unstable and to form cylindrical micelles. The φ 1 -dependence of the different membrane properties as obtained for asymmetric membranes with φ 2 = 0 is summarized in Table IV. This table contains the same numerical values as plotted in Fig. 9.

D. Comparison with previous simulation studies
Coarse-grained simulations have been previously used to study the fluid-elastic properties of bilayer membranes with two molecular components. In Ref. 37, the bilayers contained two membrane components, A and B, which consisted of a single H bead and a single chain with four and two C beads, respectively. Thus, the two membrane components differed in the length of their hydrophobic chains. Planar bilayers assembled from these lipids were studied using molecular dynamics simulations with explicit water. The bending rigidity κ was measured for different mole fractions of the binary mixture and was found to vary nonmonotonically as a function of mole fraction, qualitatively similar to the behavior found in the present study.
In Ref. 38, both membrane components consisted of a single H bead and one chain with two C beads but the H beads of the two components were different in size. The component with the large head bead had a conical shape, whereas the component with the small head bead had an inverted conical shape. These two-component membranes were assembled into spherical nanovesicles and studied by molecular dynamics simulations with implicit solvent. The curvature of the nanovesicles was varied between 1/(20 nm) and 1/(3 nm). The membrane components underwent frequent flip-flops between the two leaflets of these highly curved vesicle membranes which led, for sufficiently long run times, to certain equilibrium values of the mole fractions φ o and φ i of the large-head component in the outer and inner leaflet. These equilibrated mole fractions varied with the mean curvature 1/R of the spherical vesicles and satisfied the relations 38 which provides an example for curvature-induced sorting of the two membrane components.
By contrast, the simulations described here investigate the spontaneous curvature induced in weakly curved bilayers by a compositional asymmetry between the two leaflets as described by constant mole fractions φ 1 and φ 2 of the upper and lower leaflet. Indeed, the two lipid components SH and LH underwent essentially no flip-flops and the compositional asymmetry remained unchanged on the time scales of our simulations which had a typical run time of 20 µs. As a result, we find a linear dependence of the spontaneous curvature m on the compositional asymmetry φ 1 − φ 2 , see Eq. (39), which is very different from the logarithmic mole fraction dependence in Eq. (40) as obtained in Ref. 38.

E. Comparison with experiments on POPC/GM1 bilayers
As mentioned in the introduction, the glycolipid GM1 provides one example for a large-head lipid. The spontaneous curvature generated by GM1 has been recently studied experimentally for giant unilamellar vesicles (GUVs) formed by the phospholipid POPC with small amounts of GM1. 6,7 These GUVs had a linear size of the order of 10 µm and thus a mean curvature of the order of 1/(10 µm) which is completely negligible on nanoscopic scales. Therefore, our planar bilayers can be viewed as small segments of the weakly curved GUV membranes. In both experimental studies, the GUVs were prepared by electroformation but the details of the preparation protocol were somewhat different. Furthermore, two different methods were used to estimate the spontaneous curvature of the vesicle membranes. In Ref. 6, GUVs with spontaneously formed nanotubes were aspirated by micropipettes and the tubes were retracted by applying small suction pressures. In Ref. 7, on the other hand, larger suction pressures were used to aspirate the GUVs and, in addition, local forces were applied by optical tweezers to pull nanotubes from the GUV membranes.
The composition controlled in the tube retraction and tube pulling experiments is the overall mole fraction Φ of GM1 as used for the vesicle preparation. If all GM1 molecules were contained in the vesicle membranes, the overall mole fraction would be related to the leaflet mole fractions φ 1 and φ 2 via Φ = 1 2 (φ 1 + φ 2 ). However, a significant fraction of the GM1 molecules remains in the aqueous solution which implies the inequality between the leaflet mole fractions φ 1 and φ 2 and the overall mole fraction Φ used in the vesicle preparation.
In the tube pulling experiments, 7 the mole fractions φ 1 and φ 2 in the outer and inner membrane leaflets of the GUV membranes have been estimated via electroporation and fluorescence microscopy. As a result, the leaflet mole fractions φ 2 Φ and φ 1 Φ/4 have been obtained for overall mole fraction Φ which implies 1 2 (φ 1 + φ 2 ) 5Φ/8 < Φ, in agreement with the inequality in Eq. (41). Furthermore, the spontaneous curvature m as deduced experimentally by tube pulling was found to vary linearly with the mole fraction difference φ 1 − φ 2 , in agreement with our simulation results; see Eq. (39) and Fig. 9(c). The prefactor 1/(0.63 me ) in Eq. (39), which applies to our SH/LH bilayers, is about 2.5 times larger than the prefactor obtained from the tube pulling experiments on POPC/GM1 bilayers as studied in Ref. 7. In the tube retraction experiments, 6 the mole fractions φ 1 and φ 2 of the bilayer leaflets have not been measured but the spontaneous curvature values deduced from the latter experiments were about twice as large as those from the tube pulling method 7 for the same overall mole fraction Φ. Therefore, the spontaneous curvatures obtained for our SH/LH bilayers are quite comparable to those deduced from the tube retraction experiments on POPC/GM1 bilayers in Ref. 6.

F. Outlook on membrane buds and necks
The spontaneous curvature has a strong impact on the morphology of membranes and vesicles. One particularly interesting example is provided by the formation of buds which are connected to the mother vesicle by closed membrane necks. Within the spontaneous curvature model, the closure of the neck is described by the condition [39][40][41] This threshold for the formation of out-buds makes it difficult to study these buds in molecular simulations which are necessarily limited to relatively small vesicles with size R ve 20 nm. Indeed, for such a small vesicle, the inequality in Eq. (43) implies that the spontaneous curvature must exceed the threshold value m min 1/(14 nm), a condition that cannot always be achieved. The asymmetric adsorption layers of small solutes as studied in Ref. 9, for example, could generate curvatures only up to about 1/(24 nm).
On the other hand, the compositional asymmetry studied here can generate large spontaneous curvatures with m > 1/(14 nm) as follows from Eq. (39) for leaflet mole fractions with φ 1 − φ 2 > 0.18. Preliminary simulations confirm this conclusion; see the snapshots in Fig. 10. Because the planar membrane segment has the mean curvature M 1 0, the neck closure condition in Eq. (42) implies that the size of the bud is about R bu 1/(2m). For φ 1 = 0.5 and φ 2 = 0, the relationship in Eq. (39) leads to the spontaneous curvature m = 1/(1.3 me ) and a bud radius R bu 1.3 me . A simulation snapshot of such a bud with a closed neck is displayed in Fig. 10(a). If we consider a section through the bud parallel to the planar membrane segment, we obtain an exterior bud diameter of about three times the membrane thickness. Because the exterior diameter of the neck is about twice the membrane thickness, the bud attains an elongated non-spherical shape.
The neck closure condition in Eq. (42) and the threshold value in Eq. (43) are derived for the spontaneous curvature model. The latter model applies to bilayer membranes with at least one membrane component such as cholesterol that undergoes frequent flip-flops between the two leaflets. If flip-flops are rare, on the other hand, one has to take into account that the number of molecules is conserved within each leaflet. This constraint leads to an extension of the spontaneous curvature model which is then supplemented by an areadifference-elasticity term. 42,43 Our preliminary simulations indicate, however, that the frequency for flip-flops depends on the membrane curvature: on the time scales of the simulations, flip-flops are quite rare in weakly curved membrane segments but rather frequent close to the membrane neck which has a large negative Gaussian curvature. These flip-flops lead to the numerous large-head lipids which are visible in Fig. 10(a) within the lower leaflet.
Another process that involves the formation of membrane necks is the engulfment of nanoparticles. One example is displayed in Fig. 10(b). Initially, a rigid nanoparticle of radius 5d is assembled from nanoparticle (P) beads and put into contact with a pre-assembled symmetric bilayer. As the system relaxes toward its equilibrium state, the nanoparticle is completely engulfed by the nanoparticle. In this example, we used the parameter values f PH = 30 for the force between a nanoparticle bead and a head group bead and f PW = 100 for the force between a nanoparticle bead and a water bead. The latter parameter choice correspond to strong adhesion of the particle to the membrane. In the latter case, complete engulfment can be achieved even for symmetric bilayers with zero spontaneous curvature as follows from the general stability conditions derived in Ref. 44. These stability conditions also predict a strong dependence of the engulfment process on the values of the spontaneous curvature. The latter dependence is currently investigated by molecular simulations using the SH/LH bilayers introduced here.

VII. SUMMARY AND CONCLUSION
We studied bilayer membranes consisting of a binary mixture of lipids with small and large head groups as depicted in Fig. 1 by molecular simulations, using dissipative particle dynamics. The composition of the membranes was described by the mole fractions φ 1 and φ 2 of the large-head (LH) lipids in the upper and lower leaflet of the bilayer. Both symmetric membranes with φ 1 = φ 2 = φ le and asymmetric membranes with φ 1 φ 2 were investigated. We used the stress profiles s(z) across the membranes to identify the tensionless states for which the mechanical tension Σ as given by Eq. (11) vanishes. For these tensionless membranes, we measured the density and stress profiles (Figs. 2 and 3) as well as the composition dependence of the average molecular area (Fig. 4) for symmetric bilayers and of the leaflet tensions for asymmetric ones (Fig. 5). Furthermore, we determined the composition dependence of the area compressibility modulus K A that describes the response of symmetric and tensionless membranes to a small mechanical tension (Fig. 6). The numerical values of this elastic modulus are given in Table III and plotted in Fig. 8 as a function of the mole fraction φ le . We also measured the bilayer thickness me and determined the bending rigidity κ sy via the relationship (30) for different mole fractions (Table III and Fig. 8).
The curvature elasticity of asymmetric membranes is determined by their spontaneous curvature m and their bending rigidity κ. We varied the mole fraction φ 1 in the upper leaflet, for fixed mole fraction φ 2 = 0 in the lower leaflet, and determined (i) the parameter combination 2κm from the first moment of the stress profile as in Eq. (34) and (ii) the bending rigidity κ from the compositional dependence of the rigidity κ sy of symmetric membranes using the relation (36). Combining these two results, we obtained the spontaneous curvature m as a function of φ 1 ; see Table IV and Fig. 9. The simulation data are well fitted by the linear relation (37) which can be generalized to m = (φ 1 − φ 2 )/(2.5 nm) as follows from Eq. (39) with the membrane thickness me = 4 nm. The latter relation leads to the large spontaneous curvature m = 1/(63 nm) for the small bilayer asymmetry φ 1 − φ 2 = 0.04 and to the huge spontaneous curvature m = 1/(13 nm) for φ 1 − φ 2 = 0.2.
As explained in Subsection VI F, these large m-values are useful in order to study membrane processes such as bud formation and nanoparticle engulfment by molecular simulations. Two simulation snapshots of SH/LH bilayers that illustrate these two membrane processes are displayed in Fig. 10. Important issues that can be addressed by such simulations are the curvature-dependent flip-flops as visible in Fig. 10(a) and the dependence of the particle engulfment as shown in Fig. 10(b) on the spontaneous curvature. Another interesting aspect that is accessible to simulation studies using SH/LH bilayers is the global stability of membrane necks under mechanical perturbations.

ACKNOWLEDGMENTS
This study has been supported by the MaxSynBio consortium, jointly funded by the Max Planck Society and the Federal Ministry of Education and Research (BMBF).