Bending elasticity of vesicle membranes studied by Monte Carlo simulations of vesicle thermal shape fluctuations

Abstract

The membrane bending stiffness of nearly spherical lipid vesicles can be deduced from the analysis of their thermal shape fluctuations. The theoretical basis of this analysis [Milner and Safran, Phys. Rev. A: At., Mol., Opt. Phys., 1987, 36, 4371–4379] uses the mean field approximation. In this work we apply Monte Carlo simulations and estimate the error in the determination of the bending stiffness due to the approximations applied in the theory. It is less than 10%. The method presented in this work can also be used to determine the changes of the bending stiffness of biological membranes due to their chemical and/or structural modifications.

---

1 Introduction

Biomembranes are one of the main building blocks of living matter. Due to their important role in biological systems, there is a large interest in their structure and functioning principles.

A widely accepted description of the biomembrane structure is given by the model of Singer and Nicolson¹ representing the biomembrane as a lipid bilayer embedded with integral proteins. One of the factors determining the functioning of biomembranes are their mechanical properties, which play an important role in defining intermembrane interactions, membrane fusion, the motion of cells in flow, etc. The mechanical properties of biomembranes are determined to a great extent by the mechanical properties of their lipid bilayers. The lipid bilayer can be considered as a simplified model of the biomembrane and theoretical and experimental investigations of the mechanical properties of lipid bilayers are continuously rising.^2–5

In most cases the thickness of the lipid bilayer is much smaller than the typical dimension of the studied object (cell, vesicle, etc.) and the bilayer can be considered as a two-dimensional structure. The macroscopic mechanical properties of such membranes are determined by their stretching, bending, and shear moduli.⁶ In this paper we consider a special case of membranes with zero shear modulus, corresponding to a two-dimensional liquid. We assume that during the characteristic observation time, i.e. during the measurement, the number of amphiphilic molecules in each of the two monolayers comprising the bilayer does not change and that the volume enclosed by the membrane remains constant. This permits us to consider the studied object as an equilibrium structure. Note, however, that we do not impose a restriction for the membrane to be tensionless.

One of the most commonly used experimental methods for determining the bending stiffness of lipid membranes is the analysis of the thermal shape fluctuations of a nearly spherical lipid vesicle.^7–15 The theoretical basis of this analysis was proposed by Milner and Safran¹⁶ by considering the fluctuating vesicle as an ensemble of non-interacting harmonic oscillators. The time mean squares of the amplitudes of these oscillators depend on several factors, including the lateral stretching of the membrane that induces a lateral membrane tension. When such time and position dependent membrane tension is taken into account, the Hamiltonian of the fluctuating vesicle is not a function of an ensemble of independent oscillators. In order to obtain that, a mean field approximation is used, where the fluctuations of the membrane tension are neglected and the fluctuating tension is replaced with its mean value. A question arises, whether this approximation provides adequate results. Recently, an essay on this topic was done by Bivas and Tonchev,¹⁷ who used the Bogolyubov inequalities and the method of the approximating Hamiltonian.¹⁸ Necessary conditions assuring the validity of the Milner and Safran theory were found. One of these conditions is the high enough value of the mean lateral tension of the vesicle membrane.

One of the aims of the presented work is the verification of the results of the Milner and Safran theory by means of Monte Carlo simulations of the shape fluctuations of a nearly spherical lipid vesicle. The results of the simulations are considered as experimental data whose accuracy increases with the increase in the length of the simulation. The obtained results show that the determination of the bending stiffness is possible using the theory of Milner and Safran with acceptable precision.

The analysis of simulations presented in this work can serve as a tool to quantitatively measure the changes of the bending stiffness of biological membranes due to the chemical and/or structural modifications of lipid bilayers. This can be, for example, the change of the composition of lipid molecules in the bilayer (multicomponent lipid bilayers), the insertion of inclusions (like peptides and other proteins) in the bilayer,¹⁹ or the polymer coating of the vesicle (like PEGylated^20,21 and polyelectrolyte-grafted^22,23 vesicles).

2 The model

2.1 Bending energy

We consider thermal fluctuations of a lipid vesicle in thermodynamic equilibrium. The membrane of the vesicle is a fluid lipid bilayer. For the bending energy W_b of the membrane we use the standard Helfrich expression⁶ for a tensionless membrane with a zero spontaneous curvature and a fixed topology (the contribution of the Gaussian curvature to the bending energy does not depend on the fluctuations):where κ is the bending stiffness of the membrane, c₁ and c₂ are the principal curvatures of the vesicle membrane at the point under consideration and the integration is performed over the membrane area A.

The lipid bilayer is on a timescale of thermal fluctuations impermeable for water molecules and due to the low compressibility of water we can assume the vesicle's volume to be constant during thermal fluctuations.

With thermal fluctuations some lateral stretching of the membrane occurs on the scale of phospholipid molecules, however, the energy required to significantly change the area of the membrane greatly exceeds the thermal energy kT (product of the Boltzmann constant and the absolute temperature), therefore we can assume that the overall area A of the membrane remains almost constant during thermal fluctuations (ΔA ≪ A).

2.2 Randomly triangulated surface

The membrane of the vesicle is described by a set of N vertices linked with bonds of flexible length d to form a closed, randomly triangulated, self-avoiding network.^24–26 Lengths of the bonds can vary between their minimal value d_min and maximal value d_max. All vertices experience a hard-core repulsive potential at their mutual distances d_min.

In our simulations the initial state of triangulation is a pentagonal dipyramid with all the edges divided into equilateral bonds so that the network is composed of 3(N − 2) bonds forming 2(N − 2) triangles. The system is initially thermalized – evolved into the nearly spherical equilibrium state using the Metropolis–Hastings algorithm (as described below).

The thermal fluctuations of the vesicle membrane are simulated by employing the Monte Carlo method, where Monte Carlo steps are vertex moves, assuring shape fluctuations, and bond flips, assuring lateral fluidity within the membrane. The vertex move includes a random displacement of the vertex within a sphere with radius s positioned at the center of the vertex. In this work we choose d_max/d_min = 1.7 and s/d_min = 0.15, resulting in a self-avoiding nearly spherical network with approximately one-half of vertex moves accepted. The bond flip involves four vertices of the two neighboring triangles. The bond shared by the neighboring triangles is cut and reestablished between the other two, previously unconnected, vertices.

The microstates of the vesicle are sampled according to the Metropolis–Hastings algorithm. To obtain the canonical ensemble representing the system in a thermodynamical equilibrium, each individual Monte Carlo step (vertex move or bond flip) is accepted with probability min[1,exp(−ΔE/kT)], where ΔE is the energy change due to the vertex move or bond flip.

Some discussion was done in the past on choosing an appropriate discretization of the bending energy on a triangulated surface (for example, see Section 2 in Gompper and Kroll, 1996,²⁵ eqn (70) in Gompper and Kroll, 2004,²⁴ or Ramakrishnan et al., 2011²⁶). In this work we used for the discretization of the bending energy (eqn (1)) the relation^27,28

where the outer summation runs over all vertices and the inner summations run over all their nearest neighbors, R_i is the radial vector of vertex i, d_ij is the distance between vertices i and j,is the area of the cell in the dual lattice²⁷ in vertex i. Here σ_ij = d_ij[cot(θ₁) + cot(θ₂)]/2 is the distance between vertices in the dual lattice, θ₁ and θ₂ being opposite angles to side ij in the two triangles that share the common bond ij.

The volume of the vesicle is kept constant at the given value V₀ by the constraint |V − V₀| < ε_V, where ε_V must be small enough to fulfill the condition ε_V ≪ V₀, but still not so small to completely suppress the out-of-plane shape fluctuations of the membrane. The choice of ε_V depends on the discretization and is in our work taken to be the volume of the tetrahedron consisting of equilateral triangles with areas A₀/N_t, where A₀ is the area of the spherical vesicle with volume V₀ and N_t the number of triangles in the triangulated surface:

The evolution of the system is measured in Monte Carlo sweeps (mcs). One mcs consists of individual attempts to displace each of the N vertices, followed by 3N attempts to flip a randomly chosen bond. Let us just note that the bond flip to vertex move attempt ratio is connected to the lateral diffusion coefficient within the membrane, that is, to the membrane viscosity.^29,30 Diffusion also introduces a real timescale in the simulations and allows simulations of the dynamics of the modeled system, which is not considered in this work.

2.3 Obtaining elastic properties through spectral analysis of thermal fluctuations

Using the theory of Milner and Safran¹⁶ the bending stiffness K_c of the membrane can be measured from the spectral analysis of thermal fluctuations of our randomly triangulated surface of the nearly spherical vesicle. Note that the bending stiffness κ is an input parameter in our simulations and that we used a different symbol (K_c) for the bending stiffness obtained from the spectral analysis of thermal fluctuations. From now on, to distinguish the two values κ and K_c, we name them the input bending stiffness and the measured bending stiffness, respectively.

Let us consider our triangulated nearly spherical vesicle with volume V₀ and let R₀ be the radius of a sphere with the same volume. The length of the radial vector R_i(t) = R(ϑ_i,φ_i,t) from the origin to the vertex i at time t is then defined as

R(ϑ_i,φ_i,t) = R₀[1 + r(ϑ_i,φ_i,t)], (5)

where ϑ_i and φ_i are the spherical coordinates of the i-th vertex and r(ϑ_i,φ_i,t) is the relative displacement of the i-th vertex.

Relative displacements r(ϑ_i,φ_i,t) are decomposed into a series with respect to the spherical harmonics Y^m_l(ϑ_i,φ_i):

where cutoff l_max is of the order of R₀/d_min and the spherical harmonics are defined asusing associated Legendre polynomials P^m_l.

The complex coefficients u^m_l(t) can then be calculated using the relation

where the integration runs over the solid angle Ω of the sphere. The discretization of the above expression can be done aswhere Ω_i(t) is the solid angle corresponding to vertex i and the sum runs over all vertices of the triangulated surface.

The mean squared amplitudes of spherical harmonics 〈|u^m_l|²〉 are calculated by averaging the |u^m_l(t)|² values over an ensemble of microstates of the vesicle in the thermal equilibrium. Using the expression of Milner and Safran,¹⁶

the bending stiffness K_c and the dimensionless mean lateral tension [small sigma, Greek, macron] of the membrane can be obtained.

Note that eqn (10) were deduced for vesicles the membranes of which are lipid bilayers with a zero spontaneous curvature as well as for emulsion droplets the membranes of which are monolayers consisting of amphiphilic molecules with a zero spontaneous curvature.¹⁶ In the case of lipid bilayers K_c is the bending stiffness at free flip-flop¹⁵ (for the definition of bending elasticity at blocked and free flip-flop, see Helfrich⁶). If the stretching of the membrane of the emulsion droplet is high enough, then the sufficient condition for the validity of eqn (10) will be fulfilled.¹⁷ The same is valid for the stretched enough vesicle membrane.

Since the rhs of eqn (10) do not depend on the order of spherical harmonics m, the mean squared amplitudes of spherical harmonics obtained from simulations are first averaged over m and then the obtained values 〈|u_l|²〉 are used on the lhs of eqn (10):

To obtain the bending stiffness K_c and the dimensionless mean lateral tension [small sigma, Greek, macron] of the membrane together with their standard errors, the 〈|u_l|²〉 from simulations are fitted with the formula of Milner and Safran (eqn (11)) using an inverse squared variance weighted nonlinear fit.

3 Results and discussion

For each set of parameters the system is initially thermalized into a nearly spherical vesicle and then the volume is fixed. The squared amplitudes of spherical harmonics |u^m_l|² are obtained from Monte Carlo simulations as described in Section 2.

To obtain the ensemble of microstates that are statistically independent, the autocorrelations of squared amplitudes f(|u^m_l|²,τ) are calculated with the autocorrelation function

where the sums run over the discrete “time” t denoting consecutive microstates and T is the number of microstates used in the calculation of the mean 〈x〉. Let us define the decay time of |u^m_l|² as the value of τ when the autocorrelation function f(|u^m_l|²,τ) falls bellow 1/e.

Fig. 1 shows the autocorrelation functions of the few lowest relevant modes for a system with N = 1127 vertices and input bending stiffness κ = 20kT. It can be seen that the longest decay times are for |u^m₂|²i.e. the decay time decreases with the increasing degree of the spherical harmonics l, as expected. Let us denote the largest decay time of all the relevant modes for a given system with N vertices as τ_N (in Fig. 1 we have τ₁₁₂₇ ≈ 60 000). The largest decay time τ_N decreases with the increasing input bending stiffness κ, while it increases with the number of vertices N in the triangulated network.

Fig. 1 Autocorrelation functions f(|u^m_l|²,τ) of the lowest relevant degrees of spherical harmonics l = 2 and l = 3, for N = 1127 and κ = 20kT. The dashed gray horizontal line indicates the value 1/e. The decay time of a given mode is defined as time when the autocorrelation function falls below this value.

The decay time τ_N is important for our spectral analysis since it can be used to estimate the “time” interval between two microstates that can be regarded as statistically uncorrelated. The ensemble of statistically uncorrelated states is needed for the estimation of the standard errors together with the means of squared amplitudes of spherical harmonics. Those standard errors of 〈|u^m_l|²〉 have to be taken into account in the fitting procedure in eqn (11), to obtain relevant values of the bending stiffness K_c and the dimensionless mean lateral tension [small sigma, Greek, macron] of the membrane. In this work the interval between two consecutive microstates in an ensemble of statistically uncorrelated states, i.e. between two consecutive “measurements”, was always larger than three times the largest decay time τ_N. In Fig. 1, for example, the x-axis spans the “time” interval between consecutive measurements for a system with N = 1127 and κ = 20kT.

When squared amplitudes |u^m_l|² are averaged over the ensemble of microstates, the obtained 〈|u^m_l|²〉 with the same order m converge towards the same value, as shown in Fig. 2 for |u^m₂|². This is in accordance with the theory of Milner and Safran (rhs of eqn (10) are independent of m). Also Fig. 3 shows that the mean squared amplitudes obtained from simulations are independent of m. Note that our previously reported³¹ inability to observe this independence of 〈|u^m_l|²〉 on m was a result of numerical errors.

Fig. 2 Mean squared amplitudes 〈|u⁰₂|²〉 (full), 〈|u¹₂|²〉 (dashed) and 〈|u²₂|²〉 (dotted) as a function of the number of statistically independent measurements used in the averaging. The input bending stiffness κ = 20kT and the membrane is triangulated with N = 1127 vertices.

Fig. 3 Mean squared amplitudes 〈|u^m_l|²〉 for l = 2, 3 and 4, obtained from 1000 measurement for a vesicle with input bending stiffness κ = 20kT and triangulated with N = 1127 vertices. The error bars indicate the standard error (standard deviation divided by the square-root of the number of measurements). Lines connect the points with the same degree l and are for the guide-of-eye only.

The measured bending stiffness K_c and the dimensionless mean lateral tension [small sigma, Greek, macron] are obtained from the mean squared amplitudes as described in Section 2.3. The result of a fitting procedure for K_c is shown in Fig. 4 as a function of the maximal degree l of spherical harmonics used in the fitting of 〈|u_l|²〉 in eqn (11) (eqn for all values from l = 2 up to the maximal degree l are taken into account).

Fig. 4 Measured bending stiffness K_c together with standard error (error bars) as a function of the maximal degree l of spherical harmonics used in the calculation of K_c and [small sigma, Greek, macron]. The value of the input bending stiffness κ = 20kT is indicated with a horizontal dashed line. The membrane is triangulated with N = 3127 vertices and 200 statistically independent microstates are measured (between each measured microstate is an interval of 2 × 10⁶ mcs). Lines connecting the points are for the guide-of-eye only.

The measured bending stiffness K_c is shown in Fig. 5 as a function of the number of vertices N of the triangulated surface. As expected, the difference between the measured bending stiffness K_c and the input bending stiffness κ = 20kT decreases as we increase the number of vertices in the triangulation (i.e. as we increase the resolution of the discretization).

Fig. 5 Measured bending stiffness K_c together with standard error (error bars) as a function of the number of vertices N used in the triangulation of the membrane, for the input bending stiffness κ = 20kT. The dimensionless mean lateral tension [small sigma, Greek, macron] for the same sets of measurements is also shown. Lines connecting the points are for the guide-of-eye only.

Fig. 5 also shows the obtained values of the dimensionless mean lateral tension [small sigma, Greek, macron] for the same sets of measurements. Let us note that the measured K_c should not depend on the value of [small sigma, Greek, macron]. The mean lateral tension in the membrane depends on the value of the fixed volume of the vesicle, i.e. how much the vesicle is “swollen”. This is somewhat arbitrarily chosen by picking a random microstate in the thermodynamical equilibrium when fixing the volume and starting the measuring procedure for K_c and [small sigma, Greek, macron]. As expected, the exact choice of the equilibrium microstate used when fixing the volume, i.e. the value of [small sigma, Greek, macron], does not observably influence the measured K_c.

Fig. 6 shows the relative difference between the measured and the input bending stiffness, K_c/κ − 1, for different values of the input bending stiffness κ. It can be seen that increasing the input bending stiffness decreases the mismatch between the input and the measured bending stiffness. Note that, as already reported above, the correlation times of squared amplitudes decrease with the increasing bending stiffness of the membrane.

Fig. 6 Relative difference between the measured and the input bending stiffness as a function of the input bending stiffness for the membrane triangulated with N = 1447 vertices.

A group of methods for the investigation of lipid–water systems exists, the results of which are explained by the fluctuations of the membranes of these systems. These are the different kinds of scattering: neutron,^32,33 laser light,³⁴etc. Zilman and Granek,^32,33 investigating the neutron scattering in the lamellar lipid–water phase, used periodic boundary conditions and developed the out-of-plane fluctuations of the layers in the Fourier series. The periodic boundary conditions resemble monodisperse vesicular suspensions, with their repeating distance being the analog of the radius of the vesicle. The authors determined the decay rate of different fluctuation modes from their experimental data. Due to the absence of real timescale in our Monte Carlo simulations, we cannot determine the true values of these decay rates from simulations and a direct comparison with the experimental results cannot be made. As already noted in Section 2.2, real timescale could be introduced in our simulations through the lateral diffusion of vertices on the randomly triangulated surface, which is beyond the scope of this work.

Brocca, et al.³⁴ have investigated the shape fluctuations of large unilamellar lipid vesicles using laser light scattering. The analysis of their experimental results permits the simultaneous deduction of the time mean squares 〈|u^m_l(t)|²〉 (see eqn (10)) of the amplitudes u^m_l(t) and the decay rate of these amplitudes. According to the theory of Milner and Safran¹⁶ a relation between these two quantities exists, however the results obtained by Brocca, et al.³⁴ do not satisfy it. This discrepancy could be explained by the fact that in the systems studied by Brocca, et al.³⁴ the requirements, assuring the validity of the theory of Milner and Safran are not fulfilled: the thickness of the membrane needs to be much less than the radius of the vesicle; the viscosity of the liquid environment inside and outside the vesicle needs to be constant; etc. It must be noted that these authors used the theory of Milner and Safran that does not take into account the dissipation of the energy due to the friction between the monolayers of the bilayer, which yields a double-exponential decay of each fluctuation mode¹⁵ instead of the mono-exponential one in the theory of Milner and Safran.¹⁶ The results of our Monte Carlo simulations cannot be used for the explanation of the discrepancy found by Brocca, et al.³⁴ As in the case when the data of Zilman and Granek^32,33 were considered, the reason is the fact that our data from the simulations do not permit the determination of the true decay rate of the fluctuation modes.

4 Conclusion

The analysis of thermally induced shape fluctuations of a nearly spherical giant lipid vesicle is one of the commonly used methods to determine the bending stiffness of the vesicle's membrane. The theoretical basis of this analysis, proposed by Milner and Safran,¹⁶ uses the mean field approximation. In the present work, the error of the determination of the bending stiffness due to the approximations used in the theory was estimated.

Monte Carlo simulations of the fluctuating nearly spherical vesicle have been performed using a randomly triangulated surface. The results for the time mean squares of the amplitudes of the fluctuations, obtained from the simulations, can be determined with an arbitrarily high precision, depending only on the length of the simulation. One of the parameters in the simulations is the input value of the bending stiffness. The obtained time mean squares of the amplitudes of the fluctuations can be considered as experimental values, which can be used for the determination of the output value of the bending stiffness by means of the theory of Milner and Safran. The theory would be “exact” if the output value of the bending stiffness would have been equal to the input one. Our results show that the difference between the two values of the bending stiffness decreases with the increase of the resolution of the triangulated network and can be well below 10%.

Therefore, we can conclude that the theory of Milner and Safran can be successfully used in the determination of the bending stiffness of the membrane of a nearly spherical lipid vesicle. According to our results, the errors due to the approximations adopted in the theory are less than 10%.

In conclusion, the analysis presented in this work can be a useful tool to predict the change of the bending stiffness of biological membranes due to their chemical modification. Altering the properties of the triangulated surface and/or introducing other membrane-interacting objects in the simulations, and then measure the change of the bending stiffness, offers many useful applications. Multicomponent lipid bilayers, membranes decorated with inclusions like peptides, polymer coated vesicles like PEGylated and polyelectrolyte-grafted vesicles, just to name a few.

Acknowledgements

The authors wish to acknowledge the Slovenian Research Agency grants (Grant No. P2-0232). MF thanks Sylvio May and Dan Kroll for stimulating discussions.

References

1. S. J. Singer and G. L. Nicolson, Science, 1972, 175, 720–731.
2. U. Seifert, Adv. Phys., 1997, 46, 13–137.
3. S. A. Safran, Adv. Phys., 1999, 48, 395–448.
4. J. F. Nagle, Faraday Discuss., 2013, 161, 11–29.
5. R. Dimova, Adv. Colloid Interface Sci., 2014, 208, 225–234.
6. W. Helfrich, Z. Naturforsch., 1973, 28, 693–703.
7. M. Schneider, J. Jenkins and W. Webb, J. Phys., 1984, 45, 1457–1472.
8. H. Engelhardt, H. Duwe and E. Sackmann, J. Phys. Lett., 1985, 46, 395–400.
9. I. Bivas, P. Hanusse, P. Bothorel, J. Lalanne and O. Aguerre-Chariol, J. Phys., 1987, 48, 855–867.
10. J. Faucon, M. D. Mitov, P. Méléard, I. Bivas and P. Bothorel, J. Phys., 1989, 50, 2389–2414.
11. T. Pott and P. Méléard, Europhys. Lett., 2002, 59, 87.
12. J. Pécréaux, H.-G. Döbereiner, J. Prost, J.-F. Joanny and P. Bassereau, Eur. Phys. J. E: Soft Matter Biol. Phys., 2004, 13, 277–290.
13. R. Rodríguez-García, L. R. Arriaga, M. Mell, L. H. Moleiro, I. López-Montero and F. Monroy, Phys. Rev. Lett., 2009, 102, 128101.
14. R. S. Gracià, N. Bezlyepkina, R. L. Knorr, R. Lipowsky and R. Dimova, Soft Matter, 2010, 6, 1472–1482.
15. I. Bivas, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 061911.
16. S. T. Milner and S. A. Safran, Phys. Rev. A: At., Mol., Opt. Phys., 1987, 36, 4371–4379.
17. I. Bivas and N. S. Tonchev, J. Phys.: Conf. Ser., 2014, 558, 012020.
18. J. G. Brankov, D. M. Danchev and N. S. Tonchev, Theory of Critical Phenomena in Finite-Size Systems: Scaling and Quantum Effects, Series in Modern Condensed Matter Physics, World Scientific, Singapore, 2000, ch. 2, vol. 9.
19. V. Kralj-Iglič, V. Heinrich, S. Svetina and B. Žekš, Eur. Phys. J. B, 1999, 10, 5–8.
20. M. Wang and M. Thanou, Pharmacol. Res., 2010, 62, 90–99.
21. D. Needham and D. H. Kim, Colloids Surf., B, 2000, 18, 183–195.
22. F. Quemeneur, M. Rinaudo, G. Maret and B. Pépin-Donat, Soft Matter, 2010, 6, 4471–4481.
23. S. Sennato, D. Truzzolillo and F. Bordi, Soft Matter, 2012, 8, 9384–9395.
24. G. Gompper and D. M. Kroll, Statistical Mechanics of Membranes and Surfaces, World Scientific, Singapore, 2004, pp. 359–426.
25. G. Gompper and D. M. Kroll, J. Phys. I, 1996, 6, 1305–1320.
26. N. Ramakrishnan, P. B. S. Kumar and J. H. Ipsen, Macromol. Theory Simul., 2011, 20, 446–450.
27. C. Itzykson, Proceedings of the GIFT Seminar, Jaca 85, World scientific, Singapore, 1986, pp. 130–188.
28. D. Espriu, Phys. Lett. B, 1987, 194, 271–276.
29. H. Noguchi and G. Gompper, Phys. Rev. Lett., 2004, 93, 258102.
30. H. Noguchi and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 011901.
31. M. Fošnarič, S. Penič, A. Iglič and I. Bivas, Advances in Planar Lipid Bilayers and Liposomes, Academic Press, 2013, vol. 17, pp. 331–357.
32. A. G. Zilman and R. Granek, Phys. Rev. Lett., 1996, 77, 4788.
33. A. G. Zilman and R. Granek, Chem. Phys., 2002, 284, 195–204.
34. P. Brocca, L. Cantù, M. Corti, E. Del Favero and S. Motta, Langmuir, 2004, 20, 2141–2148.

---