Home

Search

Collections

Journals

About

Contact us

My IOPscience

Monte Carlo simulations of fluid vesicles

This content has been downloaded from IOPscience. Please scroll down to see the full text. 2015 J. Phys.: Condens. Matter 27 273104 (http://iopscience.iop.org/0953-8984/27/27/273104) View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 157.182.150.22 This content was downloaded on 25/06/2015 at 02:10

Please note that terms and conditions apply.

Journal of Physics: Condensed Matter J. Phys.: Condens. Matter 27 (2015) 273104 (14pp)

doi:10.1088/0953-8984/27/27/273104

Topical Review

Monte Carlo simulations of fluid vesicles K K Sreeja1 , John H Ipsen2 and P B Sunil Kumar1,2 1

Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India MEMPHYS-Center for Biomembrane Physics, Department of Physics and Chemistry, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark

2

E-mail: [email protected] Received 31 January 2015, revised 24 March 2015 Accepted for publication 7 April 2015 Published 18 June 2015 Abstract

Lipid vesicles are closed two dimensional fluid surfaces that are studied extensively as model systems for understanding the physical properties of biological membranes. Here we review the recent developments in the Monte Carlo techniques for simulating fluid vesicles and discuss some of their applications. The technique, which treats the membrane as an elastic sheet, is most suitable for the study of large scale conformations of membranes. The model can be used to study vesicles with fixed and varying topologies. Here we focus on the case of multi-component membranes with the local lipid and protein composition coupled to the membrane curvature leading to a variety of shapes. The phase diagram is more intriguing in the case of fluid vesicles having an in-plane orientational order that induce anisotropic directional curvatures. Methods to explore the steady state morphological structures due to active flux of materials have also been described in the context of Monte Carlo simulations. Keywords: lipid membranes, vesicle shapes, membrane inclusions, ordered membranes, Monte Carlo simulations (Some figures may appear in colour only in the online journal)

The primary forces responsible for the formation and stability of membranes are hydrophobic interactions. Owing to these interactions amphiphiles such as phospholipids spontaneously aggregate in water to form bilayer membranes. The selfassembled nature of membranes implies that it is the bending rigidity and not the surface tension which determines its large scale shape deformations, which is key to a multitude of cellular processes including budding and fusion of transport vesicles and protein sorting. Shape deformations of membranes cover many length and time scales and can vary from simple, symmetric, primitive spherical shapes, commonly seen in the case of the plasma membrane, to the highly complex, convoluted structures displayed by organelles such as the endoplasmic reticulum and the golgi. Though the molecular players involved in determining these shapes have been extensively investigated in cell biology, a generic framework to explain all the observed shapes does not exist. It is in this context that simple physical models of membranes, which capture the large scale generic properties, become important [3, 4].

1. Introduction

Biological membranes are approximately 5 nm thick quasitwo-dimensional fluid sheets composed of a variety of lipids, proteins and sugars. They separate the inner and outer environments of the cell and compartmentalize the inner organelles of the cell. These membranes support and to some extent control an amazingly specialized protein-based machinery which is crucial for a variety of physiological functions such as transmembrane transport, cell signaling and cell motility [1]. Recent developments in biology seem to indicate that the fluid mosaic model for cell membranes, proposed by Singer and Nicolson [2], with the lipid bilayer functioning only as a medium to support the protein machinery, may be too simple to be realistic. The current understanding is that the membrane composition plays an active role in controlling the protein behavior which in turn determines the local shape and composition of the membrane. A membrane is composed of a variety of lipids, which vary in size and in the nature of their hydrocarbon chains and polar head groups. 0953-8984/15/273104+14$33.00

1

© 2015 IOP Publishing Ltd Printed in the UK

J. Phys.: Condens. Matter 27 (2015) 273104

Topical Review

The internal organelles of eukaryotic cells are subjected to a steady flux of materials(lipids/proteins) in the form of membrane bound vesicles which fuse into and fission off from them. It is important to understand how the morphology of organelles emerge as a consequence of the molecular processes and physical forces involved in this transport [5– 7]. One possibility is that the shapes of these organelles, which are highly conserved across species, are in a dynamical steady state. It is also proposed that the interaction between membranes and the cytoskeletal network can also lead to stability of convoluted shapes. Apart from the out of plane deformations of the membrane, the in-plane organization of lipids and cholesterol, known as lipid rafts [8], also play a significant role in cellular transport and cell signaling [9]. However, elucidation of the fundamental issues including the mechanisms leading to the formation of rafts, their stability and finite size remain elusive. Considering the complexity of the problem, the best way to understand the role played by different components is to build models and check the conformity of the results between experiments and simulations. It is noted that, in computer simulations it is necessary to use a model that can access the large length and time scale regime relevant to this problem. The physical properties of biological membranes, which could be easily prepared in the laboratory, are believed to be similar to those of these lipid membranes. This has triggered a lot of experimental observations on the shapes and in-plane morphology of vesicular membranes, the giant unilamellar vesicles (GUV) [10]. Such model membrane studies have been especially useful to understand the in-plane organization of lipids and the induced shape changes in multi-component lipid membranes [11–19]. Such studies have helped in testing one of the main tenets of the raft hypothesis, that these membrane domains are formed solely by lipid–lipid interactions. It is now believed that equilibrium processes arising from lipid–lipid interactions alone may not be able to stabilize the nanometer sized domains seen in biological membranes. Significant progress is also made in making complexes of lipid membranes with cytoskeletal proteins—reconstituting membrane fractions extracted from cells, with energy regenerating components. The ultimate aim of these studies is to understand the factors stabilizing complex organelle morphologies in a cell. The advantages of such a bottom-up approach is obviously that we know the effect of every component included, allowing for systematic modeling of the processes, unlike processes in biological membranes wherein the membrane composition variation itself is not easy to determine. The fact that lipid membranes are self-assembled fluid interfaces separating two bulk fluid regimes presents a major challenge in its modeling. Bending energy of the membrane, which is the major component in determining its structure, is of the order of a few tens of kB T . This means entropy plays an important role in determining the conformations of these macromolecular assemblies at all length scales [20], with both in-plane compositional fluctuations and out-of-plane shape fluctuations playing important roles. Lipid membranes are amenable to topology changes, though the energy barrier for this change is expected to be large. The bilayer nature

of the lipid membranes also plays an important role in certain processes, for example in determining the membrane mediated interaction between inclusions and in-plane lipid compositional fluctuations that modify the bilayer thickness. When it comes to understanding the dynamics of these systems, since membranes are embedded in a fluid medium, the solvent mediated hydrodynamic interactions cannot be ignored. Models for lipid membranes used today include many, but not all, features of the membrane listed above. In this review, we will focus only on the Monte Carlo models for lipid membranes. A brief comparison with other models will be presented at the end. 2. Dynamical triangulation Monte Carlo simulations

The phenomenological model proposed by Canham and Helfrich [21], has been used for decades to explore the energetics of fluid membranes. The model neglects the flip flop motion of lipids between the two layers and assumes the membrane to be a quasi two- dimensional elastic sheet, with the bending elastic energy related to the local curvature in the following form.
  κ Hmean = (1) (2H − C0 )2 + κG G dS, s 2 where H and G represent the mean curvature and Gaussian curvature. The spontaneous curvature, C0 , is the preferred curvature set by the lipids at equilibrium, which is equal to zero for symmetric membrane and can be positive or negative otherwise. The molecular details of the membrane, like lipid–lipid interaction and lipid–solvent interactions are included in the coefficients, κ the bending rigidity and κG the Gaussian rigidity. Equation (1) is the bending energy of a thin elastic sheet in the incompressible limit, where the area is assumed to be a constant, with vanishing shear modulus due the unrestricted fluidity. A more general form for the energy of the membrane, which includes the effect of surface tension, area difference between and inner and outer leaves of the membrane and the osmotic pressure difference between inside and outside of closed surface, is given by,
  κ Hsurface = (2H − C0 )2 + κG G + γ dS s 2
KA π 2 + − A + P dV . (2) (A ) 0 2 AD 2 v Here γ is the surface tension which is non-zero in a stretched or compressed membrane and vanishes in the case of self-assembled vesicles. A is the area difference between the two monolayers of the membrane and the term containing this contribution in the Hamiltonian is referred to as the area difference elasticity term [22]. A0 is the equilibrium area difference and D is the separation between the monolayers. The last term in the above expression represents the contribution to the energy resulting from the asymmetry in the chemical environment between the inside and outside of the vesicle, through the osmotic pressure difference P . The above Canham–Helfrich theory describes a bilayer with a uniform lateral distribution of lipids and proteins and 2

J. Phys.: Condens. Matter 27 (2015) 273104

Topical Review

Figure 1. Monte Carlo moves performed in DTMC with an acceptance probability Pacc . (a) Vertex move changes the vertex position form {X} to {X
} and (b) Link flip changes the triangulation from {T } to {T
}.

focuses on the midsurface of a bilayer. In reality the mechanical properties of the membrane depend on its heterogeneous chemical composition. This can be incorporated by making the model parameters depend on local composition and is discussed later in this article. However, to make connection to biological membranes one needs to determine exact local composition and the dependence of membrane parameters on it, which is not easily accessible experimentally. Dynamically triangulated surfaces have been extensively used as models for fluid membranes and vesicles wherein a discrete version of the Hamiltonian given in equation (2) is used to carry out Monte Carlo simulations. In such dynamically triangulated Monte Carlo simulations(DTMCs) the continuum surface of spherical topology is replaced by a triangulated mesh of NT = 2(Nv − 2) triangles, with Nv vertices connected by 3(Nv − 2) tethers [23]. A similar relation between Nv , NT and the number of tethers can be arrived at for surfaces of arbitrary topology. The self-avoidance of the surface is ensured by assigning a hard core spherical bead of diameter √ a to each vertex and restricting the maximal tether length to 3a. This is in general not sufficient to impose strict self-avoidance, hence a constraint on the largest angle between the normals of the two faces sharing a tether is imposed. As shown in figure 1, the DTMC procedure implements the shape relaxation and fluidity through two independent Monte Carlo moves. The first is called the vertex move, in which a randomly chosen vertex undergoes a random displacement within a small cubic box centered at the vertex. This move allows the membrane shape to relax to an equilibrium conformation. In-plane fluidity of the membrane is ensured through bond flips, in which the common edge between two neighboring triangles is removed and a new tether is added between the previously unconnected vertices of the quadrilateral. This move ensures that in-plane displacement of the vertices is not constrained by the tether connecting it to its neighbors, thus ensuring fluidity. An excellent review of the DTMC model and its variants already

exists, see [24]. In this review we will concentrate mainly on topics that are not covered in [24]. 2.1. Discretized mean curvature

To carry out simulations using a discrete surface model one obviously needs to discretize the elastic Hamiltonian. The mean curvature term in equation (1) can be discretized in different ways. Below we describe some of the commonly used schemes to discretize the mean curvature. A popularly used method for discretization of mean curvature at any vertex on the membrane surface is one proposed by Itzykson [25] and is given by, n 1  σvv
(3) rvv
H= σv v
=1 |rvv
| where v represents a vertex on the discrete surface and v
is summed over all the vertices connected to v. rvv
is the vector pointing from the vertex v and its connected neighbor v
and σvv
= |rvv
|(cot θ1v
+ cot θ2v
)/2 is the bond length in the dual lattice. θ1v
and θ2v
are the angles opposite to the bond(vv
). σv is the area associated to the vertex v and is taken as 41 nj σvv
|rvv
| . For obtuse triangles σvv
can be negative and hence a possibility exists for the area of vertex to be negative. But in practice such cases are rare and the method has been shown to work very well for fluid membranes. Another method for discretization of the total mean curvature, integrated over the area, is based on the triangle normals and is given as [26, 27],
 1 − nα .nβ , (4) H 2 ds = v

vv


where nα and nβ are the local outward normals to the triangles α and β sharing a common bond (vv
). This form suffers from the drawback that the resulting bending modulus is geometry dependent. For example, it has been shown that, to get the 3

J. Phys.: Condens. Matter 27 (2015) 273104

Topical Review

(b)

(a)

ˆ Figure 2. (a) Shows a vertex v and its one ring neighborhood. Nˆ v is the normal at vertex v, R(e) the edge vector and Nˆ (e) is the normal vector to the edge e. (b) The signed dihedral angle (e) between faces, f 1 and f 2, sharing an edge e. Reproduced from [33] with permission of publisher John Wiley and Sons.

correct curvature energy in the continuum limit, the bending modulus for a spherical surface need to be 3/2 times that of a cylinder [28]. In the above methods the mean curvature integrated over the area of the membrane was calculated directly. Recently, Ramakrishnan et al [29], used a method to calculate the principle curvatures at any vertex of a triangulated surface [30– 32] and calculated the mean curvature at a given vertex as, H =

[c1 (v) + c2 (v)] 2

operator at a vertex v can be computed as Sv (v) =

 where area of the vertex A(v) = {f }v A(f )/3 and weight factor for an edge is W (e) = Nˆ (v).Nˆ (e). Using a Householder transformation S(v) can be converted from global coordinates to local coordinates in the tangent frame resulting in a 2 × 2 minor. The eigenvalues and eigenvectors of Sv (v) are the principal curvatures c1 (v) and c2 (v) and directions eˆ1 and eˆ2 , respectively.

(5)

2.2. Discretized Gaussian curvature

It should be noted that DTMC simulations can also be used to study topology changes in the vesicles. Such simulations have been performed to study the bicontinuous structure of micro emulsions and sponge phases [34–36]. For such simulations the second term in equation (1) is important and we need to find a form for the Gaussian curvature on the triangulated surface. As is well known, the Gauss Bonnet theorem ensures that the Gaussian curvature integrated over the surface is a constant as long as the topology of the surface is unchanged and hence this term can be ignored for simulations of vesicles with fixed spherical topology. However, this term needs to be included when the Gaussian curvature modulus is a function of local composition, even when the topology of the surface is unchanged. Discretization of the Gaussian curvature turns out to be best done, even in the case wherein we calculate the two principle curvatures explicitly, as deficit angles at every vertex. The Gaussian curvature integrated over the area at every vertex v is then,
 Gv dS = 2π − θv
, (9)

(6)

where summation is over all the faces around the vertex v.

is the weight factor, chosen proportional to the area A(f ). The normal at any edge e is calculated from the face normals of the triangles f1 (e) and f2 (e), sharing edge e, as ˆ 1 (e)] + Nˆ [f2 (e)] N[f ˆ . N(e) = |Nˆ [f1 (e)] + Nˆ [f2 (e)]|

(8)

v

where the principal curvatures c1 (v) and c2 (v) are calculated from a discretized shape operator constructed in the one ring neighborhood of vertex v. We will now look at this method in more detail. An illustration of a triangulated patch and associated measures are shown in figure 2. The unit normal to the surface at the vertex v, Nˆ (v), is calculated from the face ˆ ) and face area A(f ) at face f as, normal N(f  ˆ {f } [A(f )]N (f ) ˆ N (v) =  v , | {f }v [A(f )]Nˆ (f )|

1  W (e)P(v)† Se (e)P(v), A(v) {e}

(7)

To calculate the principal curvatures at any vertex v, first a ˆ ˆ ⊗ B(e)] at an discretized shape operator Se (e) = H (e)[B(e) edge e is calculated from the signed dihedral angle H (e) =   ˆ 2|R(e)| cos( (e)/2) and the binormal B(e) = Nˆ (e) × R(e) ˆ to the edge e. Here R(e) is the unit vector along the edge e. Using a projection operator P(v) = I − Nˆ (v) ⊗ Nˆ (v), which projects Se (e) into the tangent plane at vertex v, the shape

s

4

vv


J. Phys.: Condens. Matter 27 (2015) 273104

Topical Review

Figure 3. Monte Carlo move for the in-plane field.

where v
is over all connected neighbors of a vertex and θv
is the angle at the vertex subtended by each triangle at the vertex and the sum is over all triangles with v as a common vertex. The DTMC model has been successful in reproducing the equilibrium phase diagram of the vesicles and has provided a good understanding of the role of thermal fluctuations in determining their shapes [37–42]. It was also used for understanding how the coupling between membrane shape and local compositional changes affect the membrane shape and the kinetics of phase segregation in membranes [19, 43–45]. It is known that significant deformation of the membrane shape and topology changes on vesicles can be induced by curvature inducing protein inclusions. Such inclusions are generally anisotropic in shape and can cause anisotropy in the bending modulus and preferred local curvature of the membrane. We will now discuss a modified form of the Helfrich Hamiltonian that will allow us to include in-plane orientational order and directional spontaneous curvature.

local spontaneous curvatures and κ and κ⊥ , respectively, are the bending stiffness along nˆ and nˆ ⊥ . Equation (10) can be discretized, to carry out DTMC simulations, by calculating the principal curvature at every vertex as described in section 2 [29]. An in-plane field nˆ can be defined in a local Darboux frame, given by the principal directions eˆ1 and eˆ2 and the vertex normal Nˆ (v), with ϕ being the angle between nˆ and the maximally curved direction eˆ1 . The directional curvatures on the membrane, parallel and perpendicular to the orientation n, ˆ can then be calculated using Euler’s formula as Hn, = c1 cos2 ϕ(v) + c2 sin2 ϕ(v) and Hn,⊥ = c2 cos2 ϕ(v) + c1 sin2 ϕ(v). The membrane mediated interaction between the in-plane vectors n, ˆ at any two neighboring vertices, is governed by the Lebwohl–Lasher model, which is the discrete version of the one constant Frank free energy,   3 1 cos2 φ(v, v
) − . (11) HLL = − LL 2 2

3. Beyond Helfrich model: membranes with in-plane orientational fields

In [29], the angle φ(v, v
) between any two vectors on the tangent plane at neighboring vertices v and v
is computed using a parallel transport operation. The nematic–nematic orientational interaction can be replaced by an XY-like interaction [51] when the inclusions have a polar character [33]. It is also possible to include a long-range interaction between the inclusions, directly through the embedding space. DTMC should now also take into account the orientational degree of freedom of the in-plane vector. Towards this a vertex is chosen randomly and the orientation of the vector nˆ at that vertex is rotated to a new, randomly chosen direction in the tangent plane, as shown in figure 3. This rotation of nˆ , which allows for the relaxation of the orientational order, is performed keeping the vertex positions and bonds between the vertices fixed. All moves, in this dynamic triangulation Monte Carlo with in-plane orientational order (DTMCO) are now accepted through a Metropolis scheme with the total Hamiltonian which is a sum of the mean and directional curvature contributions to the elastic energy and the in-plane orientational energy of the field n, ˆ (12) Htot = Hmean + HLL + Hanis .

Proteins and other heterogeneous membranes inclusions play an important role in determining the morphologies of functionally important membrane conformations [46]. For example, interaction of proteins with the membrane is a major determinant in stabilizing highly curved, non-symmetric organelles such as the Golgi, the endoplasmic reticulum and the inner tubes of mitochondria. Molecular simulations at the level of single proteins [47–49] have shown that the curvature induced by most membrane associated proteins have an anisotropic deformation profile, which in turn calls for the inclusion of a protein orientation dependent elastic modulus into the current model for biological membranes. A possible extension to equation (1), to include an in-plane orientational field, was proposed by Frank and Kardar [50]. In this model the anisotropic nature of inclusions is introduced as an in-plane nematic field nˆ defined on the tangent plane of the membrane. The coupling of this field with the membrane curvature is anisotropic in nature and has the form
  κ κ⊥  Hanis = (Hn, − C0 )2 + (Hn,⊥ − C0⊥ )2 dA. 2 s 2 (10)

3.1. Conformations of vesicles with nematic order

The most significant change in the vesicle conformation coupled with nematic ordering, is the absence of the branched polymer phase when the bending rigidity κ is zero. Figure 4

Hn, is the curvature along direction nˆ and Hn,⊥ is the direc tional curvature along nˆ ⊥ . C0 and C0⊥ are the corresponding 5

J. Phys.: Condens. Matter 27 (2015) 273104

Topical Review

(a)

(b) 

Figure 4. Equilibrium configuration of a vesicle decorated with in-plane field, with κ = 0, c0 = 0, LL = 3.0, κ⊥ = 0 for (a) κ = 0 and

(b) κ = 20. Reprinted figure with permission from [29] ©2010 by the American Physical Society.

(a) Oblate

(b) Tubular

(c) Disc

(d) Inner tubes

(e) Tubular

(f) Corkscrew

(g) Branched

(h) Branched



Figure 5. The shapes generated by a vesicle decorated by in plane fields. C0 = 0 (a), −0.3 (b), −0.4 (c), −0.6 (d), 0.2 (e), 0.4 (f), 0.5 (g)

and 0.6 (h) with κ = 10 kB T , κ = 5 kB T , κ⊥ = 0 and εLL = 3.0. Reproduced from [52] with permission of The Royal Society of Chemistry. 

shows the equilibrium configurations of nematic membranes with κ = 0 for two different values of the directional curvature modulus κ . It can be seen that, on the scale of the vesicle size, nematic order stiffens the membrane, wiping out the branched polymer phase. Unlike the influence of mean curvature mediated membrane stiffness the equilibrium shape is not spherical, but is a twisted pillow with a half defect at each of the corners. Inclusion of a directional rigidity further elongates the vesicle, bringing the pair of +1/2 disclinations closer. Figure 5 shows the conformations of a vesicle with mean curvature rigidity κ = 10 kB T , κ⊥ = 0 and κ = 5 kB T for  different values of C0 . It can be seen that the local directional  deformation, determined by C0 , stabilizes a variety of shapes  relevant to biological cells. For negative values of C0 , the   observed shapes include oblate (C0 = 0), tubular(C0 = −0.3),



disc (C0 = −0.4) and inner tubes (C0 = −0.6). For positive  values of C0 the vesicle transforms from a tubular shape  (C0 = 0.3) to corkscrew, a tube spiralled about its long axis   (C0 = 0.4) and then to branched (C0 > 0.5). The tube diameter and the average nematic orientation depends on the  value of C0 . 3.2. Role of defects in determining the shape of vesicles with in-plane order

It is well known that in-plane orientational order on a two-dimensional membrane surface can result in singular points where the orientational order cannot be defined. The average orientation of the field rotates by nπ along a closed loop encircling these singular points. Such topological 6

J. Phys.: Condens. Matter 27 (2015) 273104

Topical Review

0.05 0.04

κ = 5 ; C0 = 0.4

κ = 0 ; C0 = 0

κ = 5 ; C0 = 0.3

κ = 5 ; C0 = 0

P (ξ)

κ = 5 ; C0 = 0.2

+1/2

0.03 Defects at equal distance

0.02

−1/2

(a)

0

(b)

+1/2 pair

0.01 0

+1/2

separation between defects pairs.

25

50

75

100

Geodesic distance, ξ

125

pair of −1/2 Annihilation of +1/2 and −1/2

Figure 6. P (ξ ) is the distribution of the geodesic distance between

defects, for shapes corresponding to different values of C0 . Reproduced from [52] with permission of The Royal Society of Chemistry.

(c)

(d)

Figure 7. The formation of tubes and branches are driven by proliferation of a +1/2 and −1/2 defect pairs. (a) +1/2 defects promotes positively curved regions, (b) −1/2 defects favour negatively curved regions, (c) pairs of −1/2 defects stabilize branches and (d) isolated +1/2 and −1/2 defects annihilate each other. Reproduced from [52] with permission of The Royal Society of Chemistry.

singularities, which cannot be removed by local rotation of the in-plane vector, are known as disclinations [53] of strength n. In the presence of a coupling between the inplane nematic order and the membrane curvature, these defects can significantly alter the morphology of the membrane. It is known that positive defects favor positive Gaussian curvature and defects with negative strength favor negative Gaussian curvature [54]. Morphology of the surface in turn alters the interaction between the defects [55]. DTMCO simulations provide a suitable platform to study this interplay between defects and membrane morphology. DTMCO simulations of vesicles has shown that the disclinations can trigger the formation of tubes and other biologically relevant shapes [52]. Defects are formed on the vesicle surface decorated by in-plane fields due to geometric frustration. On a surface with spherical topology the total strength of the defects should be +2. Since the energy of these defects scales as n2 [56], in the absence of any coupling between orientational order and curvature of the membrane, the lowest energy configuration, on a surface of spherical topology, is four +1/2 disclinations. Since similarly charged defects repel [56] they should be placed as far away as possible. In figure 6 we show the distribution of geodesic distance ξ between defects, where a single peak(circles) represents tetrahedral shapes indicating four equidistant disclinations. The presence of anisotropic curvature alters the interaction between the defects. When κ = 0 the peak splits into two (stars), indicating two +1/2 defects coming closer to form defect pairs and the pairs moving away from each other as the vesicle deforms into an oblate shape (see figure 5(a)). As shown in figures 5(e) and (f), in the presence of a directional  spontaneous curvature C0 > 0 the defects within a pair come closer at the tip of the tubes and tube length increases. This results in the ξ peaks of P (ξ ) moving away from each other.

between the sign of the defect and the sign of the curvature, these defects of opposite signs sometimes repel each other separating the negatively curved necks from the positively curved tips of the tubes. A sequence of snapshots depicting such defect induced branching is shown in figure 7, which is obtained by quenching an equilibrated vesicle from LL = 0 to  LL = 3, keeping κ = 5, κ⊥ = 0, κ = 0 and C0 = 0.4. The necks and tubes are found to be stable when a pair of −1/2 defects are at the neck and a pair of 1/2 defects are at the tip of the tube. It is found that isolated +1/2 and −1/2 defects move towards each other and annihilate. The size of the necks  and diameter of the tubes are determined by C0 . 3.4. Aggregation of anisotropic inclusions and associated shape changes in vesicles

In the previous section we looked at complete coverage of the vesicle surface by nematic field, forcing the presence of disclinations. Significant membrane shape deformations can result even when the surface is only partially decorated by in-plane fields [57], study of which can provide insight into the mechanism of shape generation due to curvature generating molecules such as proteins with BAR domain and caveolin [58, 59]. Below we will refer to these molecules as inclusions or nematogens. Such surface deformations are mainly due to the cooperative effect arising from the ability of inclusions to curve the membrane and use this curvature to mediate aggregation on the surface. Aggregation can also result from direct interaction between inclusions, which could be orientation dependent or isotropic. The nematic orientational interaction given in equation (11) is modified here to incorporate an additional isotropic interaction between the

3.3. Thermal fluctuations and generation of defects and tubes

Thermal fluctuations can generate defect pairs of opposite signs. Due to the deformability of the surface and the coupling 7

J. Phys.: Condens. Matter 27 (2015) 273104

Topical Review



Figure 8. Equilibrium configurations of vesicle for varying concentration of nematogens with κ = 10, κ = 5, κ⊥ = 0, C0 = 0.5, C0⊥ = 0,

J = 0 and LL = 3.0. Reprinted from Publication [57] ©2013 with permission from Elsevier. (a) NF = 0.1 N. (b) NF = 0.2 N. (c) NF = 0.3 N. (d) NF = 0.4 N. (e) NF = 0.5 N. (f) NF = 0.6 N. (g) NF = 0.7 N.



Figure 9. Equilibrium configurations of vesicle for varying concentration of nematogens with κ = 10, κ = 5, κ⊥ = 0, C0 = −0.5,

C0⊥

= 0, J = 0 and LL = 3.0. Reprinted from Publication [57] ©2013 with permission from Elsevier. (a) NF = 0.1. (b) NF = 0.2. (c) NF = 0.3. (d) NF = 0.4. (e) NF = 0.5. (f) NF = 0.6. (g) NF = 0.7.

inclusions.    J 3 1 2
HLL = − LL cos φ(v, v ) − I v Iv
, 2 2 2

to 0.7 the aggregates induce complex shape transformations in the vesicle. Aggregates form nematically ordered domains with well-defined curvature characteristics such as ridges or cylindrical rims. As the concentration increases larger rims are formed which stabilize to disc like structures, as seen from figure 8 for NF = 0.1 to 0.4. At intermediate concentrations discs and tubes coexist (see NF = 0.5), and for further higher concentrations discs become unstable and tube dominated structures are seen. The effect of inclusion concentration, when the directional spontaneous curvature is negative, is shown in figure 9. The interesting phenomena observed here are tubular invaginations like the ones seen in endosomes and T-tubules of muscle cells. It has been experimentally shown that, Exo 70, a membrane binding protein which can induce negative spontaneous curvature, induces tubular membrane invaginations through an

(13) where Iv and Iv
are density variables, equal to 1 when a vertex is occupied by a nematogen and zero otherwise. The aggregation of these curvature active inclusions is different from the lipid domain formation due to the orientational interaction and the anisotropy in the curvature they induce. The inclusions can generate either a positive curvature, such as F-BAR domain proteins or a negative curvature, such as I-BAR domain proteins. The conformational changes, for a positive value of directional  spontaneous curvature C0 , with varying concentration are shown in figure 8. For a surface coverage in the range NF = 0.1 8

J. Phys.: Condens. Matter 27 (2015) 273104

Topical Review



Figure 10. Equilibrium membrane conformations (spherical, disc and tubular) with NF = 0.3, κ = 20, κ = 5, C0 = 0.5 and LL = 3 for

different ranges of J . Reprinted from Publication [57] ©2013 with permission from Elsevier. (a) Spherical. (b) Disc. (c) Tubular.

model described earlier. A more detailed analysis of this problem will be published elsewhere. The model for protein induced clustering contains two types of lipids, labelled A and B and one type of protein. A vertex variable takes value φv = ±1 depending on whether that vertex is occupied by A or B type lipids. Protein occupancy is labeled by sv , with sv = 1 when the vertex is occupied by a protein and zero otherwise. We discuss here the case wherein 30% of the lipids (i.e. 30% of vertices in the triangulated membrane) are of type A and the rest of the vertices are of B type and 20% of the vertices are occupied by proteins. The proposed model allows for the co-existence of the different types of lipids and proteins on the same vertex. The discrete Hamiltonian for protein–lipid interaction is given by,  sv φ v
. (14) Hnφ = −jnφ

oligomerization based mechanism [49]. Even at low protein concentrations (NF = 0.1–0.3), tubes growing to the interior are seen. For NF > 0.3 tubes disappear and more saddle-like surfaces appear. For larger concentrations (NF > 0.8) both saddles and tubes coexist. At low inclusion densities the orientational interaction between neighboring inclusions alone can still lead to aggregation resulting in nematic ordered domains. However, in such cases the preferred curvature and low interfacial tension results in elongated domains which forms lip-like structures, as shown in figure 10(b). Obviously, adding an isotropic repulsive interaction, J < 0, destroys this ordering and domain structure leading to an isotropic and homogeneous distribution of inclusions, see figure 10(a). When J > 0, the added interfacial tension causes the boundary of the domain to shrink and the lip deforms into a tube, as shown in figure 10(c). This effect of J in the formation of tubular structures is discussed in detail in [57].

v,v


When jnφ > 0 proteins would like to be in the vicinity of A type lipids. In the summation v represents all the vertices and v
runs over all vertices directly connected to v and the vertex v itself, which means a protein placed at a vertex interacts with lipid at that vertex and that of its one ring neighborhood. Monte Carlo procedure includes a Kawasaki exchange between A and B to allow the diffusive motion of lipids. The lipid–lipid interaction energy is of the form Hφ = −jφ φv φv
, where v
sum does not include the vertex v. The net Hamiltonian is thus the sum of these terms and the curvature energy, given in equation (12). The MC time evolution of domains and the associated shape deformations are shown in figure 11. Where there are only lipid–lipid interactions (panel (a)) the well known aggregation and budding of a two component vesicle is seen [64]. For the conformations shown in panels (b)–(d) lipid–lipid interactions are absent, i.e. jφ = 0 and only lipid– protein interaction jnφ > 0 is present. Domains rich in one type of lipid can also be seen in all these cases. In panel (b),  where jφ = LL = c0 = 0 and jnφ = 1, vesicles with only isotropic interaction between lipids and protein are present. Since there is no real phase separation, the boundaries of the domains are not well defined as in the case of (a). The budding of component A is not observed in these cases as the interfacial tension, which is the main driving force for the budding of domains, is absent in the interactions.

4. Kinetics of phase separation and shape changes in multicomponent vesicles

In this section we will extend the DTMC membrane model to include a coupling between an inclusion and the local lipid composition. The lateral organization of lipids and proteins in the bilayer is significant in cell signaling, membrane trafficking and other cellular functions. Biological membranes are known to form functionally and compositionally distinct domains whose size varies from a nanometer to micrometers. For example, the well known lipid rafts are understood to be a dynamic domain of lipids and cholesterol with specific composition [60]. The ability of lipid rafts to selectively host specific proteins while excluding others makes them suitable to organize cellular activities. The larger assemblies of lipids are formed mainly due to lipid–lipid, lipid-protein and protein–protein interactions. However, how the domains maintain their nanoscale structure and how lipids and proteins are recruited into them is far from clear. On the modeling side, the formation of lipid domains and three-dimensional budding of domains in vesicles, based on lipid–lipid interactions, has been studied previously using DTMC and other simulation methods [27, 61, 62, 63, 64]. In this section we will discuss a special case of the protein induced lipid clustering due to lipid– protein and protein–protein interactions using the DTMCO 9

J. Phys.: Condens. Matter 27 (2015) 273104

Topical Review

(a)

Ncluster

n%=0; j φ=1

(b)

n%=20; jnφ=1;

=0; c0= 0.0

n%=20; jnφ=1 ;

=1; c0= 0.0

n%=20; jnφ=1;

=1; c0= 0.6

100

(c) 10000

MC steps (t)

1e+05

Figure 12. Cluster formation kinetics in multicomponent vesicles. Average number of clusters(Ncluster ) as a function of monte carlo steps is shown. Circles represent lipid clusters with pure lipid–lipid interactions with jφ = 1 and no proteins present on it, i.e. n% = 0. Squares, diamonds and triangles show the cluster size without lipid–lipid interaction, i.e., jφ = 0 and with protein–lipid affinity jnφ = 1. In the case of diamonds and triangles in addition to protein–lipid interaction, there is an in-plane interaction between proteins, LL = 1.0. Triangles show the cluster sizes in the presence of proteins with induced curvature along only one principal direction.

(d)

Initial cluster growth with and without in-plane interaction of proteins follows approximately the same power law as observed in the case of vesicle with lipid–lipid interaction. However, the directional curvature of the proteins gives rise to a different mechanism of growth; the system displays enhanced cluster formation dynamics wherein the mean cluster size approaches its saturation value in a short time.

MC steps Figure 11. Lipid segregation and shape deformations with lipid–lipid interaction and lipid–protein interaction as a function of Monte Carlo steps. Vertices of type A lipids are shown with red spheres and all other vertices representB type lipids. Proteins are shown with black lines. (a) Vertices with zero protein coverage, i.e. n% = 0 and the clusters are formed due to lipid–lipid interactions, with jφ = 1. (b)–(d) contains 20% vertices with proteins. (b) Shows the vesicle conformations for jφ = LL = c0 = 0 and jnφ = 1, which includes explicit protein–lipid interactions. (c) Here protein–protein orientational interaction is present with jφ = c0 = 0 and jnφ = LL = 1. (d) Shows the shapes with protein induced curvature c0 = 0.6 when jφ = 0 and jnφ = LL = 1.

5. Equilibrium shapes of charged vesicles

Monte Carlo simulations can also be used to determine the configurations of charged membranes. Electrostatic interactions are key to many of the protein induced morphological changes in membranes [65, 66]. Negatively charged lipids like PS, PI(Phosphatidylinositol) and polyvalent PI species like PI(4)P, PI(4,5)P2, PI(3,4,5)P3 provide a local electrostatic potential on the membrane [67]. The oppositely charged domains in the proteins adsorb on to this surface resulting in partial or complete charge neutralization. In this section we show how DTMC techniques can be used to model the shape deformations due to the electrostatics mediated adsorption of curvature inducing proteins or macro ions on a charged vesicle system. In the model system a fraction of vertices on the discretized vesicle are assumed to carry a charge −q. An equal number of oppositely charged ions called counter ions(CIs), in the space outside the vesicle, in a periodic box of side L, ensures charge neutrality. Hard core repulsion between vertices and CIs ensures impenetrability of the CIs through the membrane. The Monte Carlo steps now also include the CI moves. CIs move diffusively in the simulation box and may condense on approaching a charged site on the membrane. It should be noted that fluidity on the membrane surface ensures the

In panel (c), in addition to Hnφ we introduce an in-plane orientational interaction to the proteins, i.e. LL = 1. The lipids segregate and the proteins get oriented, but it does not seem to affect the shapes of the vesicle. When the  proteins induce a directional curvature, i.e. when c0 = 0.6, the lipids aggregate and the vesicle tubulates, as shown in panel (d). In the previous simulations of nematic membranes the protein cluster sizes computed for a value of nematic– nematic interaction strength LL = 1 were smaller than those observed in panel (c). So it is clear from above that both inplane orientational interactions and protein-lipid interactions are responsible for the formation of clusters. Figure 12 shows the kinetics of domain coarsening. Different symbols corresponds to panels (a)–(d) in figure 11. Average cluster size with only lipid–lipid interaction is compared with t 1/3 the scaling proposed in the previous studies for the case of a two component vesicle with equal number of A and B vertices. Cluster formation is observed to be faster in cases where there are non-zero lipid–protein interactions. 10

J. Phys.: Condens. Matter 27 (2015) 273104

(a)

Topical Review

(b)

Figure 13. (a) Initial configuration of a vesicle with 50% charged vertices. (b) Equilibrium conformations of a charged vesicle as a function of electrostatic constant A with condensation induced curvature C0 = 0.3. Vertices which are represented in red color are charged but uncondensed, while vertices in green color denote charged vertices with condensed CIs. CIs are not shown for clarity.

diffusivity of charges on the membrane surface. A CI and a vertex on the membrane are marked as condensed pair when the CI distance to the charged vertex is below a fixed distance = 1.1a. The ability of a bound protein to curve the membrane is modeled through the introduction of spontaneous curvature C0 at CI condensed vertices. The parameter C0 thus couples the elastic energy to electrostatic energy. The total Hamiltonian is the sum of elastic and electrostatic contributions, where the later is the sum over Coulomb interaction between all pairs of charges in the system, given as
 qi qj κ H = {2H − C0 }2 dS + A , (15) 2 s rij i,j

% of condensed ions

80 60 40 c0=0.0 c0=0.1 c0=0.3 c0=0.5

20 0 0

where A = 4π 10 r is the electrostatic constant. r is the dielectric constant of the medium. qi and qj are the charges present at the vertices and CIs. In general Ewald summation is to be employed to evaluate the electrostatic component of the Hamiltonian [68]. However, in the regime where the CI condensation is high, a periodic minimum image calculation is found to be equally good. Equilibrium configurations of vesicles obtained from Monte Carlo simulations are given in figure 13. The shapes of the vesicle are the result of competing elastic, electrostatic and entropic contributions to the energy. In all cases the simulations are carried out from an initially quasi-spherical vesicle configuration of radius ∼12a, see figure 13(a). At low values of A(< 1.0), the condensation of CI is low and the electrostatic repulsion between the nodes keeps the vesicle in a quasi spherical shape. As the strength of electrostatic interaction increases the vesicle shape transforms from an elongated tubular structure to a branched structure. These shape transformations are induced by the CI condensation through the spontaneous curvature and the change in the effective electrostatic interaction between vertices that it induces. Figure 14 shows the percentage of ions condensed as a function of A for different values of C0 . As expected, the number of condensed counter ions increases with A. Interestingly, the amount of condensation also depends on the spontaneous curvature C0 induced by the CI. As condensation number increases total mean curvature of the vesicle increases, prompting a shape transformation to tubular and branched

2

4

A

6

8

10

Figure 14. Percentage of condensed CIs as a function of electrostatic constant A for spontaneous curvature C0 = 0.0, 0.1, 0.3, 0.5. Condensation number increases with both A and C0 .

shapes. It is clear from figure 14 that positively curved regions are more favorable for CI condensation. This is expected since the charges on the positively curved regions of the membrane are more accessible to CIs. Since CIs induce positive curvature, condensation is also a cooperative phenomenon. 6. Membrane interaction with polymers and nano particles

Apart from the stiff anisotropic inclusions discussed in the previous section, the effect of interaction between semiflexible polymers and nano particles on the membrane morphologies have also been studied using Monte Carlo simulations. Simulations of fluid vesicles that contain a single worm-like polymer chain, with varying persistence length, showed a transition of the polymer from an isotropic disordered random conformation to an ordered toroidal coil [69]. Concurrent to this polymer shape transition the vesicle adopts an oblate shape. These studies have relevance in the context of understanding the effects of cytoplasmic biopolymers on the shape of cells. Fosnaric and coworkers have looked at the wrapping of a charged colloid by a triangulated fluid membrane with fraction of its vertices being oppositely charged [70]. They 11

J. Phys.: Condens. Matter 27 (2015) 273104

Topical Review

observed a discontinuous wrapping transition from partial to complete or almost complete wrapping as a function of the electrostatics screening length and fraction of charged vertices on the membrane. Increasing the electrostatic screening length and decreasing the bending rigidity was found to weaken the transition. Interestingly, the preferred shape of the vesicle, on wrapping of the spherical colloid, becomes nonaxisymmetric when the spontaneous curvature of the charged vertices is of opposite sign to the curvature of the colloid. Another interesting problem that is studied extensively using MC simulations is the interactions between membranes and colloids. The effect of colloids on the conformation of membranes and the membrane conformation mediated assembly of colloids on its surface were investigated [70–73]. The interaction of nano particles with a membrane and how this leads to its assembly and internalization of particles are important questions for our understanding of nano-toxicity. It has been demonstrated, using DTMC simulations, that strong elastic pair interactions between pairs of particles emerge as the result of their binding to membranes. Depending on the excess area available, these interactions can lead to partially wrapped pairs to linear aggregates of particles confined to membrane tubes [72, 73]. They indicate tube formation as the dominant mechanism for internalization of nano particles [72]. Matthews and Likos looked at assembly of patchy particles on fluid membranes [74]. They explored the phase diagram of these composites by varying interaction between colloidal particles, that between colloids and membrane and the stiffness of the membrane. The interactions between patchy colloids are chosen such that the structures formed on the membrane resemble intermediate viral capsomers (referred to as cores) or clathrin coats. The final structure depends on the number of patches on the colloids and the maximum allowed angle between patchy particles before the attractive interaction between the patches decays [74]. In all cases the attractive interaction between the membrane and colloids was found to increase the probability of assembly for parameter sets wherein bulk assembly is not preferred. However, membrane bending modulus has a varying effect on the formation of cores and clathrin-like coats. For cores, membrane rigidity shows an interesting nonmonotonic dependence, being less favored for very deformable membranes and disappearing for the stiffest. On the other hand, for clathrin like particles, the promotion of assembly persists for less deformable membranes.

the membrane bound compartments could be influenced by the active out-of-equilibrium processes of fission and fusion of material. One of the main effect of these transport events is to create and remove curvature locally on the membrane with the help of a variety of enzymes that convert chemical energy to mechanical work. An obvious question is then, what is the steady state shapes of fluid membranes subjected to such non-equilibrium curvature generation and removal events? In DTMC simulations these active events of curvature generation and removal were represented by a scalar field φ at every vertex i, which takes values +1 or −1. φ in turn sets the spontaneous curvature C0 , in equation (1), at the vertex i to C0i = M0 (1 + φi )/2. φi = ±1, implies a curvature generation and removal process, respectively. Here M0 is a parameter that sets the value of preferred curvature [78]. The transition probabilities for φi
−φi are taken to be independent of each other. The explicit form of these transition rates are,  N+ 1 (16) P+→− = − N 1 + exp(ζ [N+ − N− − A0 ]) and, P−→+ = +



N− N



1 . η + exp(−ζ [N+ − N− − A0 ])

(17)

Here, + = − = denote the mean attempt rate for these non-equilibrium curvature changes. This form for transition rates ensures that the instantaneous number of vertices N± , with φi = ±1 (with N = N+ + N− , the total number of vertices), does not deviate significantly from a desired value N±0 . Note that these transition rates are entirely dependent on the preferred asymmetry parameter, A0 ≡ N+0 − N−0 and the parameter ζ which sets the scale of fluctuations in N+ . N+0 and N−0 are the steady state mean

valuesof N+ and N− ; to ensure N±

reaches N±0 we set η = 2 NN−+ − 1 in equation (17). Note that the above transition probabilities represent non-equilibrium processes and do not obey detailed balance and are included in addition to the equilibrium Metropolis moves described earlier. In addition, such curvature changes resulting from the binding and unbinding of curvature generating complexes could be cooperative and can be accounted for by an Ising–Hamiltonian, 1 J φi φj 2 i=1 j ∈

N

Hφ = −

(18)

i

A sample of the steady state shapes that result from the activity is shown in figure 15. As can be seen from the figure, even a small fraction of active sites can change the shape of the membrane significantly. The non-equilibrium steady-state nature of these shapes was further confirmed by the the fact that on secession of activity the membrane goes back to its equilibrium spherical shape [78].

7. Active membranes: modeling non-equilibrium transport on membrane surfaces

The membrane bound organelles of eukaryotic cells are distinguished by their unique morphology and chemical composition. They emerge in the face of a steady flux of membrane bound vesicles which fuse into and fission off from them. Thus these organelles, in the trafficking pathways, are dynamic membranous structures [75]. Since many studies have indicated that the time scales of these material fluxes are at least comparable to membrane relaxation times in the highly viscous environment of the cell [76, 77], which are in the order of tens of seconds, the large-scale morphology of

8. Summary, conclusions and outlook

Dynamical triangulation Monte Carlo simulations has been used extensively in the past to study finite temperature shape phase diagrams of single- and multi-component membranes. They have also been used to study driven transport of vesicles 12

J. Phys.: Condens. Matter 27 (2015) 273104

Topical Review

processes are discussed here within the framework of DTMC simulations, to explore the steady-state shapes of a membrane with non-equilibrium curvature fluctuations. Methods to estimate the principal curvatures on triangulated surfaces require further refinement to get the Gaussian curvature at any vertex more accurately. Furthermore, computationally viable methods to incorporate equilibrium and nonequilibrium area fluctuations of the membrane should be developed in order to model membranes subjected to transport of lipids and proteins. These are some of the problems that future work using DTMC will try to address. Acknowledgment

We thank N Ramakrishnan for his critical reading and constructive comments on the manuscript. We acknowledge the use of HPCE facilities at IIT Madras. Sreeja K K thanks UGC, Government of India, for financial support. References [1] Alberts B, Bray D, Lewis J, Raff M, Roberts K and Watson J D 1994 Molecular Biology of the Cell (New York: Garland) [2] Singer S J and Nicolson G L 1972 Science 175 720–31 [3] Lipowsky R and Sackmann E 1995 Structure and Dynamics of Membranes (Amsterdam: Elsevier) [4] Ramakrishnan N, Kumar P B S and Radhakrishnan R 2014 Phys. Rep. 543 1 [5] Rafelski S M and Marshall W F 2008 Nature 9 593–602 [6] Marshall W F 2011 BMC Biol. 9 57 [7] Mart´ınez-Men´arguez J A 2013 ISRN Cell Biol. 2013 1–15 [8] Simons K and Ikonen E 1997 Nature 387 569 [9] Gousset K, Wolkers W F, Tsvetkova N M, Oliver A E, Field C L, Walker N, Crowe J H and Tablin F 2002 J. Cell. Physiol. 190 117 [10] Kas J and Sackmann E 1991 Biophys. J. 60 825–44 [11] Silvius J R, del Giudice D and Lafleur M 1996 Biochemistry 35 15198 [12] Korlach J, Schwille P, Webb W W and Feigenson G W 1999 Proc. Natl Acad. Sci. 96 8461 [13] Bagatolli L A and Gratton E 2000 Biophys. J. 78 290 [14] Dietrich C, Bagatolli L A, Volovyk Z N, Thompson N L, Levi M, Jacobson K and Gratton E 2001 Biophys. J. 80 1417 [15] Veatch S L and Keller S L 2002 Phys. Rev. Lett. 89 268101 [16] Veatch S L and Keller S L 2005 Phys. Rev. Lett. 94 148101 [17] Baugart T, Das S, Webb W W and Jenkins J T 2005 Biophys. J. 89 1067 [18] Li L, Liang X, Lin M and Yang Y 2005 J. Am. Chem. Soc. 127 17996 [19] Bagatolli L A and Kumar P B S 2009 Soft Matter 5 3234 [20] Mouritsen O G 2005 Life: As a Matter of Fat (Berlin: Springer) [21] Helfrich W 1973 Z. Nat.forsch. C 28 693 [22] Miao L, Seifert U, Wortis M and D¨obereiner H G 1994 Phys. Rev. E 49 5389–407 [23] Ho J S and Baumgartner A 1989 Phys. Rev. Lett. 63 1324–4 [24] Gompper G and Kroll D 2003 Chapter 12: Triangulated-surface models of fluctuating membranes Statistical Mechanics of Membranes and Surfaces ed D Nelson et al (Singapore: World Scientific) [25] Itzykson C 1986 Proc. of the GIFT Seminar, Jaca 85 ed J Abad et al pp 130–88 [26] Kantor Y and Nelson D R 1987 Phys. Rev. Lett. 58 2774–7 [27] Kumar P B S and Rao M 1998 Phys. Rev. Lett. 80 2489–92 [28] Gompper G 1996 J. Phys. I 6 1305–20

Figure 15. Steady-state shapes of an active membrane, (a) for

= 0.1 N/MCS and J = 0, as a function of curvature-activity coupling M0 , (b) for J = 0 and M0 = 0.8, as a function of activity rate, . The side of the stomatocyte that is curved-in, is colored differently, for clarity. (c)Steady-state shapes at = 0.1 N/MCS and M0 = 0.8, as a function of cooperativity J between active species. All configurations are obtained with κ = 20 and N+0 = 0.1 N. The locations of the active protein complexes are shown by the shaded regions.

through narrow pores and budding induced by crystalline patches on the membrane. Excellent reviews on these topics already exist [24, 43]. Instead, in this article we have reviewed some of the recent developments and applications of DTMC to understand the morphology of biological membranes. Anisotropic bending elasticity [50] induced due to inplane ordering and the presence of inclusions, is an interesting new dimension added to DTMC simulations. This model can also be used to study the interplay between orientational ordering of inclusions and the lipid composition, which may have relevance to lipid rafts in biological membranes. It should be noted that both curvature generation and curvature sensing are important properties of protein membrane interactions. It is shown that curvature sensing leads to the aggregation of inclusions, which in turn stabilizes morphologies that are otherwise untenable. Our preliminary results show that the sensitivity of membrane curvature to the binding of proteins could also be due to the interplay between electrostatic interactions and curvature generation. One of the main interests in studying membrane morphology is to understand the organization of cellular organelles. Many of these organelle membranes are subjected to a steady flux of transport vesicles and proteins, which fuse and detach from their surface through active processes. These 13

J. Phys.: Condens. Matter 27 (2015) 273104

Topical Review

[29] Ramakrishnan N, Sunil Kumar P B and Ipsen J H 2010 Phys. Rev. E 81 041922 [30] Polthier K 2002 Polyhedral surfaces of constant mean curvature PhD Thesis TU-Berlin [31] Taubin G 1995 Estimating the tensor of curvature of a surface from a polyhedral approximation Proc. Fifth Int. Conf. on Computer Vision, 1995 pp 902–7 [32] Hildebrandt K and Polthier K 2004 Eurographics 23 391–400 [33] Ramakrishnan N, Kumar P B S and Ipsen J H 2011 Macromol. Theory Simul. 20 446–50 [34] Gompper G and Kroll D M 1998 Phys. Rev. Lett. 81 2284–7 [35] Jeppesen C and Ipsen J H 1993 Eur. Phys. Lett. 22 713 [36] Peltomaeki M, Gompper G and Kroll D M 2012 J. Chem. Phys. 136 134708 [37] Gompper G and Kroll D M 1994 Phys. Rev. Lett. 73 2139–42 [38] Kroll D M and Gompper G 1995 Phys. Rev. E 51 514 [39] Kroll D M and Gompper G 1992 Science 255 968–71 [40] Boal D and Rao M 1992 Phys. Rev. A 45 R6947–50 [41] Noguchi H and Gompper G 2004 Phys. Rev. Lett. 93 258102 [42] Noguchi H and Gompper G 2005 Proc. Natl Acad. Sci. USA 102 14159–64 [43] Ramachandran S, Laradji M and Kumar P B S 2009 J. Phys. Soc. Japan 78 041006 [44] Kohyama T and Gompper G 2007 Phys. Rev. Lett. 98 198101 [45] Penic S, Mesarec L, Fosnaric M, Iglic V K, Kralj S, Gozdz W and Iglic A 2014 J. Phys.: Conf. Ser. 558 012010 [46] Atilgan E and Sun S X 2007 J. Chem. Phys. 126 095102 [47] Blood P D and Voth G A 2006 Proc. Natl Acad. Sci. 103 15068–72 (Preprint) [48] Yin Y, Arkhipov A and Schulten K 2009 Structure 17 882–92 [49] Zhao Y et al 2013 Dev. Cell 26 266–78 [50] Frank J R and Kardar M 2008 Phys. Rev. E 77 041705 [51] Chaikin P M and Lubensky T C 2000 Principles of Condensed Matter Physics 1st edn (Cambridge: Cambridge University Press) [52] Ramakrishnan N, Ipsen J H and Kumar P B S 2012 Soft Matter 8 3058–61 [53] Lavrentovich O 2001 Encyclopedia of Materials: Science and Technology (New York: Elsevier) pp 6071–6 [54] Bowick M, Nelson D R and Travesset A 2004 Phys. Rev. E 69 041102

[55] MacKintosh F and Lubensky T 1991 Phys. Rev. Lett. 67 1169–72 [56] De Gennes P G and Prost J 1993 The Physics of Liquid Crystals (The International Series of Monographs on Physics) (Oxford: Clarendon) [57] Ramakrishnan N, Sunil Kumar P B and Ipsen J 2013 Biophys. J. 104 1018–28 [58] Frost A, Unger V M and De Camilli P 2009 Cell 137 191–6 [59] Williams T M and Lisanti M P 2004 Genome Biol. 5 214 [60] Rao M and Mayor S 2014 Curr. Opin. Cell Biol. 29 126–32 [61] Laradji M and Kumar P B S 2004 Phys. Rev. Lett. 93 198105 [62] Laradji M and Kumar P B S 2004 Phys. Rev. Lett. 93 198105 [63] Kumar P B S, Gompper G and Lipowsky R 1999 Phys. Rev. E 60 4610 [64] Sunil Kumar P B, Gompper G and Lipowsky R 2001 Phys. Rev. Lett. 86 3911–4 [65] McMahon H T and Gallop J L 2005 Nature 438 590–6 [66] Li J, Zhang H, Qiu F, Yang Y and Chen J Z Y 2015 Soft Matter 11 1788–93 [67] McLaughlin S and Murray D 2005 Nature 438 605–11 [68] Frenkel D and Smit B 2001 Understanding Molecular Simulation: From Algorithms to Applications 2nd edn (New York: Academic) [69] Fosnaric M, Iglic A, Kroll D M and May S 2013 Soft Matter 9 3976–84 [70] Fosnaric M, Iglic A, Kroll D M and May S 2009 J. Chem. Phys. 131 105103 [71] Saric A and Cacciuto A 2013 Soft Matter 9 6677–95 [72] Saric A and Cacciuto A 2012 Phys. Rev. Lett. 109 188101 [73] Bahrami A H, Lipowsky R and Weikl T R 2012 Phys. Rev. Lett. 109 188102 [74] Matthews R and Likos C N 2012 Phys. Rev. Lett. 109 178302 [75] Alberts B, Johnson A, Lewis J, Raff M, Roberts K and Walter P 1994 Molecular Biology of the Cell 3rd edn (Singapore: Garland Publishing) [76] Wirtz D 2009 Annu. Rev. Biophys. 38 301–26 [77] Wieland F T, Gleason M L, Serafini T A and Rothman J E 1987 Cell 50 289–300 [78] Ramakrishnan N, Ipsen J H, Rao M and Kumar P B S 2015 Soft Matter 11 2387–93

14