Membrane simulation models from nm to $\mu$m scale

Recent developments in lipid membrane models for simulations are reviewed. To reduce computational costs, various coarse-grained molecular models have been proposed. Among them, implicit solvent (solvent-free) molecular models are relatively more coarse-grained and efficient for simulating large bilayer membranes. On a $\mu$m scale, the molecular details are typically negligible and the membrane can be described as a continuous curved surface. The theoretical models for fluid and elastic membranes with mesh or meshless discretizations are presented. As examples of applications, the dynamics of vesicles in flows, vesicle formation, and membrane fusion are presented.


I. INTRODUCTION
An amphiphilic molecule, such as a lipid and detergent, consists of hydrophilic ('water-loving') and hydrophobic ('water-fearing') parts. In aqueous solution, these molecules self-assemble into various structures depending on the relative size of the hydrophilic part; spherical or cylindrical-like micelles, bilayer membranes, and inverted micelles. In particular, the bilayer membrane of phospholipids is the basic structure of the plasma membrane and intracellular compartments of living cells, where the membranes are in a fluid phase and lipid molecules can diffuse in two-dimensional space. A vesicle (closed membrane) is considered to be a simple model of cells and also has applications as a drug-delivery system.
Many interesting behaviors of bilayer membranes have been studied such as fluid-to-gel phase transition, lateral phase separation, protein insertion, the lysis of lipid vesicles, and membrane fusion. The thickness of a lipid membrane is 5nm and cells are ∼ 10µm in diameter. Thus, the length scale of these phenomena varies from nm to µm. To investigate the bilayer structure and molecular interactions, information on a nm scale is necessary. For example, the length mismatch of 1nm between a protein and lipid can change the stability and interactions between proteins [1]. On the other hand, the morphologies of lipid vesicles and cells on a µm scale can be understood by continuous surface theories [2,3,4].
Various types of membrane models have been proposed for the simulations. They are classified in two groups depending on whether the bilayer structure is implicitly or explicitly taken into account. In the first group [Figs. 1(d), (e)], the bilayer membrane is described as a smoothly-curved mathematical surface [2,3,4]. This assumption is valid on length scales greater than the membrane thickness of 5nm. Thus, a unit segment represents not a lipid molecule, but a membrane patch consisting of hundreds to thousands of lipid molecules. The information about the bilayer properties is only reflected in this case in the values of the elastic parameters. A mesh discretization with a triangular or square mesh is typically employed for the simulation. To allow large deformations  [16]. (c) solvent-free molecular model [81]. (d) meshless membrane model [70]. (e) dynamically-triangulated membrane model [5]. Typical lengths of unit segments are shown in parentheses. Figures (a) and (b) were provided by Shinoda. of a fluid membrane, remeshing [5,6,7] is necessary. In particular, a dynamically-triangulated membrane model [5,6] is widely used for the membrane with thermal fluctuations. An alternative model is a meshless membrane, in which particles self-assemble into a membrane by potential interactions. The meshless models are very suitable for studying membrane dynamics accompanied by topological changes.
In the second group [Figs. 1(a)-(c)], amphiphilic molecules are modeled by atomistic or coarse-grained (CG) molecular models, and the solvent is taken into account either explicitly or implicitly. Although computer technology is rapidly growing, 50 ns dynamics of hundreds of lipid molecules are typical levels of recent simulations by the all-atom models. Thus, CG models have been developed to reduce the computational costs. CG molecular models with the Lenard-Jones potential [8] FIG. 2: (Color online) Bond flip of triangulated mesh to produce membrane fluidity. and with dissipative particle dynamics (DPD) [9,10] potential [11] have been widely used for the explicit-solvent simulations. Recently, the potential parameters in the CG molecular models are tuned by atomistic simulations [12,13,14,15,16], see Figs. 1(a) and (b). To further reduce the computational costs, larger CG segments (two or three segment particles per amphiphilic molecule) are employed and the solvent is implicitly represented by an effective attractive potential between the hydrophobic segments, see Figs. 1(c). This type of implicit solvent model is typically called the solvent-free model.
There is a length-scale gap between these molecular models and the curved-surface models. For the molecular models, the size of the unit segments can be gradually varied from 1Å to 1nm. For the curved-surface models, a unit segment represents 5nm, and the thermal undulations of smaller lengths are neglected.
Recently, the CG molecular simulations and triangulated membrane were reviewed in Refs. [17,18,19] and in Ref. [5], respectively. The phase separation in multicomponent membranes was reported by Kumar in this JSPSJ special issue. In this paper, we review the recent developments in the solvent-free molecular models and curved surface models [Figs. 1(c)-(e)] for singlecomponent membranes and also describes some new results. In Secs. II and III, the triangulated and meshless membrane models are described, respectively. The areadifference elasticity (ADE) model is newly applied to the dynamically-triangulated membrane. The membrane dynamics in flows are presented for both models. In Sec. IV, the solvent-free molecular models are described. The membrane fusion and the formation of polyhedral vesicles are presented as applications.

A. fluid membrane
The curvature energy of a fluid membrane is given by [2,20,21] where C 1 and C 2 are the principal curvatures at each point in the membrane. The coefficients κ andκ are the bending rigidity and saddle-splay modulus, respectively. The spontaneous curvature C 0 vanishes when lipids symmetrically distribute in both monolayers of the bilayer. Lipid membranes typically have κ = 20k B T , where k B T is the thermal energy [3]. For membranes of fixed topology without edges, the integral over the Gaussian curvature C 1 C 2 is an invariant and their properties are independent ofκ. Since the critical micelle concentration (CMC) of lipids is very low, the number of lipid molecules in a membrane is essentially constant over the typical experimental time scales. Also, the osmotic pressure generated by ions or macromolecules in solution, which cannot penetrate the lipid bilayer, keeps the internal volume essentially constant. Under the constraints of a constant volume V and constant surface area A, stomatocyte, discocyte, and prolate shapes give the energy minimum of F cv with C 0 = 0 for 0 < V * 0.59, 0.59 V * 0.65, and 0.65 V * < 1, respectively [22], where the reduced volume V * = 3V /(4πR 3 0 ) with R 0 = (A/4π) 1/2 . In a dynamically-triangulated surface model [5,6,23,24,25] of vesicles and cells, the membrane is described by N vertices which are connected by bonds (tethers) to form a triangular network. The vertices have excluded volumes. The curvature energy is discretized using angles of the neighboring faces [5,6,26] or using dual lattices of triangulation [5,6,27]. To model the fluidity of the membrane, bonds can be flipped to the diagonal of two adjacent triangles using the Monte Carlo (MC) method, see Fig. 2. The membrane viscosity η mb can be varied by the frequency of the bond flip [28,29]. This model is widely used to simulate a vesicle with spherical topology.
The membrane with open edges was also studied [24,25]. Topological changes can be taken into account by a mesh reconnection [5,30]. However, this reconnection is a discrete procedure, so that it cannot smoothly treat the dynamics.
Recently, the dynamically-triangulated model were applied to the dynamics of a fluid vesicle in flows by combining with boundary element methods [31,32] or multiparticle collision (MPC) dynamics [28,29,33,34]. MPC is a particle-based simulation method proposed by Malevanets and Kapral [35]. In a simple shear flow, v =γye x , the vesicle shows three types of motions [28,29,31,33]. (i) Tank-treading motion: The vesicle has the constant inclination angle θ with respect to the flow direction, and instead, the membrane is rotating. (ii) Tumbling: The angle θ rotates. (iii) Swinging: The angle θ oscillates around zero. At a low shear rateγ, the vesicle transits from the tank-treading to tumbling with the increasing viscosity contrast η in /η 0 or membrane viscosity η mb , where η in and η 0 are the viscosities of the internal and external fluids. At the higher shear rateγ, a swinging phase appears between the tank-treading and tumbling phases. These dynamics were also observed in experiments [36,37], and are understood by the theory for an ellipsoidal vesicle [33,38] or the perturbation theory for a quasi-spherical vesicle [39,40,41]. We found shearinduced shape transitions; elongating transitions from stomatocyte and discocyte to prolate, as well as a shrinking transition from a tumbling prolate to a tank-treading discocyte [28,29]. In capillary flow, the fluid vesicle also transits from discocyte to prolate with the increasing flow velocity [34].

B. area difference elasticity
In experiments, several morphologies of lipid vesicles were observed in the addition to the above shapes: pear, peal-neckless, and branched starfish-like shapes [3,4,42,43]. To explain them, the area difference elasticity (ADE) model is widely used [4,44]. In a vesicle, the area of two monolayers of the bilayer membrane are different with ∆A = h (C 1 + C 2 )dA, where h is the distance between two monolayers. Since the flip-flop of lipids (traverse motion between monolayers) is very slow, the area difference ∆A 0 = (N out − N in )a 0 preferred by lipids is typically different from ∆A, where N out and N in are the number of lipids in the outer and inner monolayers, respectively, and a 0 is the area per lipid. In the ADE model, the energy of this mismatch ∆A − ∆A 0 is given by with the averaged The curvature normalized by m = 4π of a spherical vesicle, is denoted by * , as m * = m/4π. We simulated vesicles using the dynamically-triangulated model for F = F cv + F ade with the harmonic constraint potentials of A and V . Figure 3 shows snapshots of the Brownian dynamics (BD) simulation, where the motions of the membrane vertices are given by the underdamped Langevin equations.
Vesicles of the ADE model also exhibit shape transitions in shear flow. A budded vesicle [see the bottom snapshots in Fig. 4] shows a tumbling motion at the low shear rateγ = 0.92 [see red lines in Fig. 4], even for η in /η 0 = 1 and η mb = 0, where prolate and oblate vesicles at any V * only show tank-treading. A similar tumbling motion was reported for a budded two-component vesicle [45]. At the higher shear rateγ = 3.68, the vesicle is found to transit into a prolate shape, which shows tank-treading [see blue lines in Fig. 4]. As the shear rate decreases, the vesicle transits back to a budded shape. The vesicle dynamics of other shapes, such as starfishes, are an interesting problem for further studies.

C. elastic membrane
The membranes (shells) of synthetic capsules, virus, and red blood cells (RBCs), have the shear elasticity µ. We call them elastic membranes. The elastic membrane can be modeled as a fixed mesh (network). The relative importance of the bending and shear elasticities is determined by the dimensionless Föppl-von Kármán number γ = ER 2 0 /κ, where E is the two-dimensional Young modulus [46]. The elastic membrane with 12 disclinations has a spherical shape at γ 150 and an icosahedral shape at γ 200. The bending or shear elasticity is dominant for the former or latter condition, respectively. The properties of viral capsids are well explained by the elastic membrane models [26,46,47].
The models of an RBC membrane have been intensively studied [34,48,49,50,51,52,53]. In the molecular-scale models [49,50], the spectorin network attached on the bilayer membrane is modeled by bonds with a worm-like chain (WLC) potential. In the continuum models, the Skalak model [48] is widely used. Recently, the comparison with the optical tweezer experiment [54] became a standard method to tune the elastic ηin/η0 = 1, η mb = 0, and N = 500. The vesicle exhibits shape transitions from a budded shape to prolate atγ = 3.68 and from prolate to the budded shape atγ = 0.92. The asphericity is the degree of deviation from a spherical shape [131] are the eigenvalues of the gyration tensor of the vesicle and R 2 g = λ1 + λ2 + λ3. The intrinsic time unit is τ = η0R 3 0 /κ. The other simulation parameters are described in Ref. [29]. parameters [50,51,52,53]. Figure 5 shows the experimental data [54] and our simulation data. In the simulation, the RBC membrane is modeled by triangular networks with 578 vertices, which are connected by a bond potential U bond = (k 1 /2)(r − r 0 ) 2 {1 + (k 2 /2)(r/r 0 − 1) 2 } with µ = ( √ 3/4)k 1 = 6 × 10 −6 N/m and κ = 2 × 10 −19 J. The nonlinear term with k 2 = 1 gives better agreement with the experimental data, but an arbitrariness in the choice of the nonlinear function still remains. The shape transitions of RBC with a stomatocyte-discocyteechinocyte sequence can be reproduced by the elastic membrane with the ADE model [55].
The dynamics of the RBCs and elastic capsules have been simulated by the elastic membrane combined with boundary element methods [56,57,58], immersed boundary methods [51,59,60,61], MPC [34], and DPD [53]. Figure 6 shows the elastic vesicles in capillary flow [34]. At the small flow velocities, the symmetry axis of the discocyte is found not to be oriented perpendicular to the cylinder axis. With the increasing flow velocity, the elastic vesicle transits into a parachute-like shape while the fluid vesicle transits into a prolate shape. The transition velocity linearly depends on the elasticities µ and κ. The results are in good agreement with the experimental data [62]. In simple shear flow, elastic capsules [61,63,64] and RBCs [53,65] transit from tumbling to tank-treading with shape oscillation (swinging) as the shear rate increases. These dynamics are caused by an energy barrier for the tank-treading rotation [66].

III. MESHLESS MEMBRANE
The meshless models are particle-based methods, which do not need bond connections between neighbor- ing particles. All of the proposed meshless methods are solvent-fee methods, but can be extended to explicit solvent versions. The first meshless model was proposed by Drouffe et al. in 1991 [67]. The particles possess an orientational degree of freedom and interact with each other via three potentials: a soft-core repulsion, an anisotropic attraction, and a hydrophobic multibody interaction. The particles self-assemble into a fluid membrane. Very recently, the modified models with pairwise potentials were pro- Recently, we proposed an alternative meshless model [70] in which particles possess no internal degrees of freedom -unlike the other models. In our model, the particles interact with each other via the potential which consists of a repulsive soft-core potential U rep with a diameter σ, an attractive potential U att , and a curvature potential U α . All three potentials only depend on the positions r i of the particles. Two types of curvature potentials [70] were proposed on the basis of the moving least-squares (MLS) method [71,72]. Here, we only describe the simpler potential, which drives the particles onto a plane: U α = i α pl (r i ). The degree of deviation from a plane, the aplanarity α pl , is defined as where λ 1 , λ 2 , and λ 3 are eigenvalues of the weighted gyration tensor, a αβ = j (α j − α G )(β j − β G )w mls (r i,j ) with α, β = x, y, z and a smoothly-truncated Gaussian function w mls (r) [70]. In this model, κ and Γ can essentially be independently varied, see Fig. 7. The bending rigidity κ linearly increases with k α and is almost independent of ε. On the other hand, the line tension Γ linearly increases with ε and is almost independent of k α for k α /k B T 10. The saddle-splay modulus is roughly estimated asκ ≃ −κ.
Although the hydrodynamic interactions are not present in the model itself, these interactions can be taken into account by combining with a particle-based hydrodynamic method such as MPC [35] and DPD [9,10]. Thus, the hydrodynamic interactions can be easily switched on or off. Static equilibrium properties can be investigated by BD or MC, which requires much less computational time. One of the disadvantages of meshless methods, compared with the mesh models, is that the vesicle volume cannot be constrained, since the exact volume is difficult to calculate. However, the volume can be naturally fixed if an explicit solvent is introduced.
The effects of the hydrodynamic interactions for the self-assembly and dissolution are investigated by comparing MPC with BD [73]. The hydrodynamic interactions are found to speed up the dynamics in both cases. The particles self-assemble into discoidal clusters, and clusters aggregate, and then large clusters form vesicles. To quantitatively investigate the closing dynamics to a vesicle, the simulations of the initially flat membranes were performed. For a membrane slightly above the critical size 4N 0 , all of them close into vesicles via a bowl-like shape. At N = N 0 , the line tension energy of a flat disk equals the bending energy of a spherical vesicle [74]. For much larger membranes of N ≫ 4N 0 , however, the membrane stochastically forms two vesicles via an S-shaped conformation, see Fig. 8. The large line tension induces a negative surface tension, and the resulting buckled shape can grow. The closing vesicle has an oblate shape in the BD simulations with large N [see the left-bottom snapshot in Fig. 8]. With the hydrodynamic interactions, the pressure of the internal fluid makes this closing shape more spherical.
When the line tension is reduced, a vesicle dissolves via an opened pore on the membrane with N (t) = N (0) − ct/2 for both MPC and BD [73]. A similar pore-opening was observed in the lysis of a lipid vesicle by detergents [75].
In simple shear flow, a spherical vesicle of the meshless membrane elongates into a prolate shape with its decreasing volume. Similar elongation is observed for the vesicles consisting of surfactants [76]. At higher shears, the vesicle ruptures into two pieces, see Fig. 9. The membrane rupture is a fundamental process for the formation of multi-lamellar vesicles in shear flows [77,78,79].

A. model
The simulations of lipid molecules with an explicit solvent require the calculation for a large number of water molecules in addition to lipid molecules. A small patch of flat membranes can be simulated with 30 water molecules per lipid by the atomistic models [80]. However, more water molecules are needed for the simulations of vesicles, since the formation of a vesicle (Fig. 8) requires large solvent space to prevent membrane interactions through the periodic boundary conditions of the simulation box. The self-assembly of amphiphilic molecules in dilute solutions requires many more water molecules. The solventfree models are more efficient tools for the simulations which require a larger solvent space. A similar solventfree approach is also frequently used in the simulations of polymers and proteins.
A common feature of the models is the requirement of an attractive potential between hydrophobic segments. This attraction mimics the "hydrophobic" interaction (hydrophobic segments dislike to contact with water). One of the simple ideas to estimate this interaction is that the hydration energy is assumed to be proportional to the solvent accessible surface area (ASA) [90,91]. Since the calculation of the ASA is a numerically time-consuming task, an effective potential, which is a function of the local density ρ of hydrophobic segments, were proposed for protein [92] and lipid molecules [81]. At a low density ρ < ρ * − 1, the potential acts as a pairwise attractive potential, but at a high density ρ > ρ * , this attraction vanishes, in which the segments are assumed to be completely surrounded by other hydrophobic segments. This multibody potential can produce a very fast lateral diffusion of molecules and a wide fluid-phase range (0.1 k B T /ε 0.9) [81]. A fluid membrane is also generated by pairwise attractive potentials, but their CMC is very high, and the membrane often coexists with isolated monomers. Although the range of the fluid phase becomes larger for wider pairwise potentials [85], their ranges are still relatively small. Figure 10 shows that  [97]. The red spheres and yellow cylinders represent the hydrophilic and hydrophobic segments of amphiphilic molecules, respectively. the lateral diffusion on the membrane becomes faster with the decreasing ρ * , where ρ * = ∞ corresponds to a pairwise potential. The attractive pairwise potential of other solvent-free models can be extended to densitydependent potentials in order to obtain faster dynamics.
The solvent-free molecular models have been applied to a variety of phenomena: self-assembly to vesicles [81,87], membrane fusion [93,94,95], membrane fission [94,96,97], the formation of polyhedral vesicles [97], pore formation [82,98], the adhesion of nanoparticles [94,99], the fluid-gel phase transition [88,100], phase separation of lipids [85], protein inclusion in membrane [98,100], and DNA-membrane complexes [101]. Figure 11 shows a polyhedral vesicle. When the bending rigidity is large compared with the vesicle radius, the lined defects at the edges of polyhedrons are induced by the mismatch of the membrane curvature with the spontaneous curvature C 0 of the monolayer. The number of edges and vertices increases with the increasing C 0 [97]. Recently, a similar defect was proposed as a kink structure in a symmetric ripple phase [102].
In early years, an alternative implicit-solvent molecular model, in which the head segment is constrained on a plane, was employed [103,104]. In this model, the membrane cannot be deformed or self-assembled. However, the diffusive motion of molecules is still present on a fixed membrane geometry and it may provide a good reference state to investigate the effects of the thermal undulations of a membrane. Recently, a phantom-solvent model was proposed [102,105], where an explicit but very simple solvent is employed. In this model, solvent particles have a repulsive interaction with amphiphilic molecules, but do not interact with each other. Thus, the solvent is in the ideal-gas equation of state, P V = nk B T , so that its pressure is easily controlled.

B. membrane fusion and fission
The membrane fusion and fission are key events in various intra-and intercellular processes, such as protein trafficking, fertilization, and viral infection. In an endocytosis pathway, a small vesicle pinches off from the plasma membrane and fuses with a lysosome. Several fusion mechanisms and fusion intermediate structures have been proposed for the fusion of biological and lipid membranes [106]. Among them, the stalk model [18,106,107,108,109,110,111] is widely accepted and is qualitatively supported by experimental studies. The first intermediate in this model, a stalk, is a hourglass-like structure that connects only the outer monolayers of the vesicles [ Fig. 12(b)]. In the stalk models, two types of pore-opening pathways from the stalk intermediate were previously proposed. Radial expansion of the stalk results in contact between the inner monolayers inside the stalk, a trans-monolayer contact state [ Fig. 12(c)]. In the original model [108], expansion of the contact area results in a disk-shaped bilayer consisting of both inner monolayers, called a hemifusion diaphragm [ Fig. 12(e)], and a fusion pore is formed in the hemifusion diaphragm [  was estimated to be too high (∼ 200k B T ) in Ref. [109], the barrier is reduced by a more quantitative estimation of the monolayer geometry [107,111]. The molecular tilt in the monolayers was proposed as the mechanism to fill the void space [110]. Recently, the free energy landscape of the fusion was also studied using the self-consistent field theory [18,112,113].
The first molecular simulation of membrane fusion was performed by a solvent-free model [93]. Small vesicles spontaneously fuse via the modified stalk pathway [ The adhesion of a nanoparticle was found to promote this fusion process [94]. Recently, the leakage between the interior and exterior of a vesicle was observed in an experiment [114]. This leakage supports the sidepore pathway. The vesicles pinched by particles form the trans-monolayer contact [ Fig. 13 [96]. On the other hand, the particle adhesion induces pore opening on the necked membrane and the pore expansion leads to the stalk formation [ Fig. 12:(f)→(d)→(b)] [94]. In large spontaneous curvature of the monolayers induces the vesicle fission via pore-opening at the lined defects [97].
Recently, membrane fusion was simulated by various molecular models such as the bond-fluctuation lattice model [115,116], explicit-solvent CG models [117,118,119,120], DPD models [121,122,123], and the atomistic model [124]. These simulations also show fusion pathways as shown in Fig. 12. The hemifusion diaphragm was found to also be formed from the side-pore of the stalk [117,119,123]. The membrane fission was also simulated in two-component vesicles [45,125,126]. Although the dependence on external forces [95] and surface tension [18,116,122,123] was investigated, conditions to determine the pathways are not fully understood. The membrane fusion and fission are induced by proteins in living cells. Their molecular mechanisms are not well understood and have not yet been simulated.

V. SUMMARY
We have presented three types of membrane models: the triangulated membrane, meshless membrane, and solvent-free molecular models. The first two models are constructed for large scale (µm) and the last model is for a molecular scale (nm). To study the bending deformation of membranes, the details on a molecular scale are typically not necessary, therefore, the mesh and meshless models are suitable and much more efficient than the molecular models. On the other hand, to study the phenomena including the non-bilayer structure, such as membrane fusion and protein insertion, the molecular scale cannot be neglected.
When compared with the phenomena in thermal equilibrium, the nonequilibrium phenomena are not well understood. In this paper, we presented several dynamic behaviors on a µm-scale in flows. On the other hand, the dynamics on a molecular scale were much less explored while a few simulation studies [127,128] were reported. The electric field is another interesting external field, which can open pores on a membrane [129] and induce shape deformation of a vesicle [130]. In living cells, the biomembrane is in a complex nonequilibrium environment. The membrane models presented in this paper are powerful tools to study the membrane behaviors under nonequilibrium as well as under equilibrium conditions.