# Numerical Study of Membrane Configurations

Authors: Luka Mesarec; Miha Fošnarič; Samo Penič; Veronika Kralj Iglič; Samo Kralj; Wojciech Góźdź; Aleš Iglič
Journal: Advances in Condensed Matter Physics
 (2014)

---

## Abstract

We studied biological membranes of spherical topology within the framework of the spontaneous curvature model. Both Monte Carlo simulations and the numerical minimization of the curvature energy were used to obtain the shapes of the vesicles. The shapes of the vesicles and their energy were calculated for different values of the reduced volume. The vesicles which exhibit in-plane ordering were also studied. Minimal models have been developed in order to study the orientational ordering in colloids coated with a thin sheet of nematic liquid crystal (nematic shells). The topological defects are always present on the surfaces with the topology of a sphere. The location of the topological defects depends strongly on the curvature of the surface. We studied the nematic ordering and the formation of topological defects on vesicles obtained by the minimization of the spontaneous curvature energy.


---

## Body

## 1. Introduction

Human red blood cells are one of the most intriguing systems in nature. The shapes of the red blood cells have been studied by many theoretical methods. The spontaneous curvature model, developed by Deuling and Helfrich [1], described the stomatocyte and discocyte shapes of the red blood cells for the first time. The shapes were calculated by solving the Euler-Lagrange equations numerically. The results of the theoretical calculations were in a very good agreement with experimental observations.The biological membrane forms a wall, which surrounds the cell and intercellular organelles [2]. It takes part in transport of nutrients, encapsulation of larger particles or viruses [3, 4], cell-to-cell communication, waste control, and many other biological processes. It is a complex system composed of lipids, carbohydrates, proteins, and many other biologically active components [5]. Lipids which form a bilayer [6] are the main building blocks of biological membranes. The membrane shape changes due to the thermal motion of the molecules which form the bilayer [7]. The thickness of the lipid bilayer is of the order of 3–5 nm. The size of the vesicles formed by the bilayer is of the order of hundreds of micrometres. Therefore, it is fully justified to use the elastic continuum approach in the theoretical description of the membrane surface.The vesicles which exhibit in-plane ordering are of particular interest. An example of such a vesicle is a colloidal particle coated with a thin sheet of nematic liquid crystal, called nematic shell [8, 9]. For such systems, minimal models were developed, which capture the main phenomena related to their orientational order. Liquid crystal molecules are oriented within the tangent plane of the shell. On nematic shells with the topology of a sphere, topological defects are always present [10]. At the origin of topological defects, orientational order is melted. The location of topological defects is strongly curvature dependent [9, 11–14]. On a spherical surface, the equilibrium configuration typically has four topological defects located in vertices of a tetrahedron in order to maximize their mutual separation, which was proved by theoretical calculations [15] and confirmed experimentally [16].The paper is organized as follows. In Section2, we introduce the theoretical models and briefly discuss the numerical procedures used to calculate the equilibrium vesicle shapes and the nematic ordering on those vesicles. In Section 3, we present the results of the numerical calculations. The vesicles calculated by two different methods are compared. The shapes of the vesicles obtained in the framework of the spontaneous curvature model are used to calculate the nematic ordering on a vesicle. The summary and conclusions are presented in Section 4.


## 2. Models and Methods

The vesicle shapes have been studied within the framework of Helfrich spontaneous curvature model [17]. We have used both Monte Carlo simulations and numerical functional minimization procedures to calculate the shapes of the vesicles. The results obtained in both methods are compared in the following section. The nematic ordering on vesicles has been studied within the two-dimensional Landau-de Gennes tensorial formalism. The equilibrium textures were calculated in Monte Carlo simulations.Many models were formulated in order to study the shape transformations of biological membranes. The model developed by Helfrich [17] is one of the most often used and the most successful. In this model, the lipid bilayer forms a homogeneous two-dimensional fluid. It is approximated by a mathematical surface and the energy of the membrane depends on the mean and Gaussian curvatures of that surface. The local bending energy density of the membrane can be written as [1, 17–20]
(1)fb=κ2C1+C2-C02+κ¯C1C2,
where κ and κ¯ are bending constants, C1 and C2 are principal curvatures of a vesicle surface, and C0 is the spontaneous curvature. Equation (1) is a limiting case of a more general model for bending energy, which also includes different membrane components and was described in [21, 22]. Overall bending energy Ftot is calculated by integrating (1) over the entire surface of a vesicle:
(2)Ftot=∫S‍fbdS,
where dS is an infinitesimal element of the vesicle area S. The Gauss-Bonnet theorem states that the last term on the right-hand side of (1), integrated over the entire surface of a vesicle, is constant for closed surfaces of fixed topology. That term only contributes a constant to the bending energy and is not taken into account in our calculations, because our surface has a fixed topology. The bending elasticity modulus κ was experimentally measured in [23] (see also [24] and references therein). Similar models were proposed by Canham [25] and Evans [18].


### 2.1. Minimization Procedure

We assumed that the vesicle surface is the surface of revolution, with rotational symmetry about thez-axis. To describe such vesicles, we only need to define a vesicle profile on a plane (Figure 1). The curve that defines the vesicle profile is then rotated around the z-axis by an angle φ=2π. The surface of the vesicle is constructed in this manner. To describe the vesicle profile on the r-z plane, we introduce the arc length s of the profile curve and an angle θ(s). θ(s) is the angle of the tangent to the profile curve with the plane that is perpendicular to the axis of rotation z (Figure 1). If the function θ(s) is known, the vesicle profile can be calculated by the following parametric equations:
(3)z(s)=∫0s‍sinθ(s′)ds′,rs=∫0s‍cos⁡θs′ds′,
where r(s) and z(s) are the coordinates of the vesicle profile in the r-z plane. We can expand the function θ(s) in the Fourier series [26]:
(4)θ(s)=θ0Lss+∑n=1N‍ansinnπLss,
where Ls is the profile length, N is the number of Fourier modes, and an are the Fourier amplitudes, which are calculated when the bending energy Ftot is minimized. For closed shapes, we apply the following boundary conditions: θ(0)=0, θ(Ls)=π, and r(0)=r(Ls)=0, which means that θ0=π in (4). The bending energy Ftot is a function of Fourier amplitudes an and the profile length Ls; therefore we need to do the numerical minimization of the function of many variables in order to calculate vesicle shapes [26–28]. In the minimization process, constraints on the surface area and the volume of the vesicle are imposed in order to set a fixed value of the reduced volume v. The reduced volume v is defined as the ratio of the vesicle volume to the volume of the sphere with the same surface area as a given vesicle.Figure 1
Vesicle profile representation onr-z plane.


### 2.2. Monte Carlo Simulations

The fluid vesicle is represented by a set ofN vertices that are linked by tethers (i.e., bonds) of variable length l so as to form a closed, randomly triangulated, self-avoiding surface [29, 30]. The lengths of the tethers can vary between a minimal value, lmin⁡, and a maximal value, lmax⁡. The self-avoidance of the surface is implemented by ensuring that no vertex can penetrate through the triangulated surface and that no bond can cut through another bond. In our simulations, we use lmax⁡/lmin⁡=1.7.The randomly triangulated network acquires its lateral fluidity from a bond flip mechanism [29, 30]. A single bond flip involves the four vertices of two neighbouring triangles. The tether connecting the two vertices in a diagonal direction is cut and reestablished between the other two, previously unconnected, vertices.The microstates of the vesicle are sampled according to the Metropolis algorithm, with the energy for a given microstate(5)E=Ftot-ΔpV,
where the first contribution is the elastic bending energy of the vesicle (see (2)) and the second contribution accounts for the energy change with the change of the volume of the vesicle, V, due to the pressure difference, Δp, inside and outside of the vesicle. Our vesicle consists of a symmetric membrane (including the absence of a mismatch between the lateral areas of the two individual membrane leaflets), so we do not need to include the spontaneous curvature (C0=0). The bending energy Ftot of the discretized vesicle (i.e., of the triangulated network) is calculated as described by Gompper and Kroll [29, 30]; for a recent review, see [31].The evolution of the system is reported inMonte Carlo sweeps (mcs). One mcs consists of individual attempts to displace each of the N vertices by a random increment in the sphere with the radius δ, centered at the vertex, followed by RbN attempts to flip a randomly chosen bond. We denote Rb as the bond-flip ratio, which defines how many attempts to flip a bond are made per one attempt to move a vertex in one mcs. Note that the bond-flip ratio is related to the lateral diffusion coefficient within the membrane, that is, to the membrane viscosity [32, 33]. The diffusion also introduces a real time scale in the simulations and allows for the simulation of the dynamics of the modelled system (not considered in this work). In our simulations, we have chosen Rb=3 and δ/lmin⁡=0.15.The investigated vesicle consists ofN=1447 vertices, which are connected with 3(N-2)=4335 bonds to form Nt=2(N-2)=2890 triangles. The vesicle, if spherical, has a radius of approximately 13. During simulations, the coordination number for each vertex (i.e., the number of its nearest neighbours, Z) is allowed to vary between 4 and 8. For the bending stiffness of the vesicle, we use κ=10kBT, where kB is the Boltzmann constant and T the absolute energy. In the following, we use lmin⁡ as the unit of length and kBT as the unit of energy.


### 2.3. Nematic Shells

We study the nematic ordering on smooth, closed, axial-symmetric surfaces, which we have calculated within the spontaneous curvature model. We use the minimal model, developed to study nematic shells [8, 34], which considers effects related to in-plane ordering within vesicles. To describe the orientational ordering of molecules, which are bound to lie on the local tangent plane on a surface, we introduce the surface order tensor Q. On the surface, we introduce a local orthonormal basis (e1,e2) in which we can represent Q as [11, 35]
(6)Q=q0(e1⊗e1-e2⊗e2)+qm(e1⊗e2+e2⊗e1),
where q0 and qm are scalar functions in the chosen coordinate system. The tensor Q can also be written in a diagonal form:
(7)Q=λn⊗n-n⊥⊗n⊥,
where λ and -λ are eigenvalues of eigenvectors n and n⊥. Value λ is a measure for the orientational order and is bound to lie on an interval λ∈[0,1/2]. The value λ=1/2 corresponds to the maximal degree of orientational order, while the value λ=0 corresponds to an isotropic state with no orientational order at all. The points on the vesicle exhibiting λ=0 usually fingerprint topological defects, since at the core of topological defects the nematic director is not defined. The key characteristic property of a topological defect is its topological charge [35, 36]. According to the theorem of Poincaré [37], the net topological charge is determined by a surface topology and it equals 2 for a sphere and all surfaces obtained by smoothly deforming a sphere, which are also the surfaces of our interest.If we know the values ofq0 and qm at a given point, we can calculate the direction of molecules and the orientational order at that point. The orientational ordering is calculated with the minimization of free energy. In the simplest form, the dimensionless free energy density is [35]
(8)f~=Rξ02tTr⁡Q~2+14Tr⁡Q~22+12∇~sQ~2,
where R is the characteristic length of a vesicle and ξ0 the nematic coherence length which is the shortest length in the system, typically in the nanoscopic range. An ordered phase can occur when the reduced temperature t is negative, which means that the system is below critical temperature. The operator ∇~s represents the surface gradient. The tilde notation is used because we have scaled the tensor Q and all the distances in order to get a dimensionless expression.Vesicles of any shape can be coated with a thin layer of nematic liquid crystal to get the previously described nematic shells. In order to calculate the total free energy, we integrate (8) over the whole surface of a vesicle. We use the Monte Carlo method to minimize the free energy. We are randomly changing values of q~0 and q~m at random points on the surface, until we reach the equilibrium configuration. The tilde notation is used, because we scaled q0 and qm in order to get a dimensionless expression for energy (see (8)).


## 2.1. Minimization Procedure

We assumed that the vesicle surface is the surface of revolution, with rotational symmetry about thez-axis. To describe such vesicles, we only need to define a vesicle profile on a plane (Figure 1). The curve that defines the vesicle profile is then rotated around the z-axis by an angle φ=2π. The surface of the vesicle is constructed in this manner. To describe the vesicle profile on the r-z plane, we introduce the arc length s of the profile curve and an angle θ(s). θ(s) is the angle of the tangent to the profile curve with the plane that is perpendicular to the axis of rotation z (Figure 1). If the function θ(s) is known, the vesicle profile can be calculated by the following parametric equations:
(3)z(s)=∫0s‍sinθ(s′)ds′,rs=∫0s‍cos⁡θs′ds′,
where r(s) and z(s) are the coordinates of the vesicle profile in the r-z plane. We can expand the function θ(s) in the Fourier series [26]:
(4)θ(s)=θ0Lss+∑n=1N‍ansinnπLss,
where Ls is the profile length, N is the number of Fourier modes, and an are the Fourier amplitudes, which are calculated when the bending energy Ftot is minimized. For closed shapes, we apply the following boundary conditions: θ(0)=0, θ(Ls)=π, and r(0)=r(Ls)=0, which means that θ0=π in (4). The bending energy Ftot is a function of Fourier amplitudes an and the profile length Ls; therefore we need to do the numerical minimization of the function of many variables in order to calculate vesicle shapes [26–28]. In the minimization process, constraints on the surface area and the volume of the vesicle are imposed in order to set a fixed value of the reduced volume v. The reduced volume v is defined as the ratio of the vesicle volume to the volume of the sphere with the same surface area as a given vesicle.Figure 1
Vesicle profile representation onr-z plane.


## 2.2. Monte Carlo Simulations

The fluid vesicle is represented by a set ofN vertices that are linked by tethers (i.e., bonds) of variable length l so as to form a closed, randomly triangulated, self-avoiding surface [29, 30]. The lengths of the tethers can vary between a minimal value, lmin⁡, and a maximal value, lmax⁡. The self-avoidance of the surface is implemented by ensuring that no vertex can penetrate through the triangulated surface and that no bond can cut through another bond. In our simulations, we use lmax⁡/lmin⁡=1.7.The randomly triangulated network acquires its lateral fluidity from a bond flip mechanism [29, 30]. A single bond flip involves the four vertices of two neighbouring triangles. The tether connecting the two vertices in a diagonal direction is cut and reestablished between the other two, previously unconnected, vertices.The microstates of the vesicle are sampled according to the Metropolis algorithm, with the energy for a given microstate(5)E=Ftot-ΔpV,
where the first contribution is the elastic bending energy of the vesicle (see (2)) and the second contribution accounts for the energy change with the change of the volume of the vesicle, V, due to the pressure difference, Δp, inside and outside of the vesicle. Our vesicle consists of a symmetric membrane (including the absence of a mismatch between the lateral areas of the two individual membrane leaflets), so we do not need to include the spontaneous curvature (C0=0). The bending energy Ftot of the discretized vesicle (i.e., of the triangulated network) is calculated as described by Gompper and Kroll [29, 30]; for a recent review, see [31].The evolution of the system is reported inMonte Carlo sweeps (mcs). One mcs consists of individual attempts to displace each of the N vertices by a random increment in the sphere with the radius δ, centered at the vertex, followed by RbN attempts to flip a randomly chosen bond. We denote Rb as the bond-flip ratio, which defines how many attempts to flip a bond are made per one attempt to move a vertex in one mcs. Note that the bond-flip ratio is related to the lateral diffusion coefficient within the membrane, that is, to the membrane viscosity [32, 33]. The diffusion also introduces a real time scale in the simulations and allows for the simulation of the dynamics of the modelled system (not considered in this work). In our simulations, we have chosen Rb=3 and δ/lmin⁡=0.15.The investigated vesicle consists ofN=1447 vertices, which are connected with 3(N-2)=4335 bonds to form Nt=2(N-2)=2890 triangles. The vesicle, if spherical, has a radius of approximately 13. During simulations, the coordination number for each vertex (i.e., the number of its nearest neighbours, Z) is allowed to vary between 4 and 8. For the bending stiffness of the vesicle, we use κ=10kBT, where kB is the Boltzmann constant and T the absolute energy. In the following, we use lmin⁡ as the unit of length and kBT as the unit of energy.


## 2.3. Nematic Shells

We study the nematic ordering on smooth, closed, axial-symmetric surfaces, which we have calculated within the spontaneous curvature model. We use the minimal model, developed to study nematic shells [8, 34], which considers effects related to in-plane ordering within vesicles. To describe the orientational ordering of molecules, which are bound to lie on the local tangent plane on a surface, we introduce the surface order tensor Q. On the surface, we introduce a local orthonormal basis (e1,e2) in which we can represent Q as [11, 35]
(6)Q=q0(e1⊗e1-e2⊗e2)+qm(e1⊗e2+e2⊗e1),
where q0 and qm are scalar functions in the chosen coordinate system. The tensor Q can also be written in a diagonal form:
(7)Q=λn⊗n-n⊥⊗n⊥,
where λ and -λ are eigenvalues of eigenvectors n and n⊥. Value λ is a measure for the orientational order and is bound to lie on an interval λ∈[0,1/2]. The value λ=1/2 corresponds to the maximal degree of orientational order, while the value λ=0 corresponds to an isotropic state with no orientational order at all. The points on the vesicle exhibiting λ=0 usually fingerprint topological defects, since at the core of topological defects the nematic director is not defined. The key characteristic property of a topological defect is its topological charge [35, 36]. According to the theorem of Poincaré [37], the net topological charge is determined by a surface topology and it equals 2 for a sphere and all surfaces obtained by smoothly deforming a sphere, which are also the surfaces of our interest.If we know the values ofq0 and qm at a given point, we can calculate the direction of molecules and the orientational order at that point. The orientational ordering is calculated with the minimization of free energy. In the simplest form, the dimensionless free energy density is [35]
(8)f~=Rξ02tTr⁡Q~2+14Tr⁡Q~22+12∇~sQ~2,
where R is the characteristic length of a vesicle and ξ0 the nematic coherence length which is the shortest length in the system, typically in the nanoscopic range. An ordered phase can occur when the reduced temperature t is negative, which means that the system is below critical temperature. The operator ∇~s represents the surface gradient. The tilde notation is used because we have scaled the tensor Q and all the distances in order to get a dimensionless expression.Vesicles of any shape can be coated with a thin layer of nematic liquid crystal to get the previously described nematic shells. In order to calculate the total free energy, we integrate (8) over the whole surface of a vesicle. We use the Monte Carlo method to minimize the free energy. We are randomly changing values of q~0 and q~m at random points on the surface, until we reach the equilibrium configuration. The tilde notation is used, because we scaled q0 and qm in order to get a dimensionless expression for energy (see (8)).


## 3. Results and Discussion

The result of minimizing the bending energyFtot, given by (2), is presented in Figure 2. Three classes of vesicle shapes in their ground states (i.e., in states with the minimal Ftot) are shown along with their energies for different values of the reduced volume v. We have obtained stomatocyte, oblate, and prolate shapes, each of them being the equilibrium shapes on some range of the reduced volume. Similar results were obtained in [38, 39].Figure 2
Bending energy of closed vesicles in unitsF=Ftot/(8πκ) as a function of reduced volume v. Shapes were calculated within the spontaneous curvature model for C0=0.In Figure3, we compare the shapes obtained with two different methods described in Sections 2.1 and 2.2. The method described in Section 2.1 allows us to set the value for the reduced volume v, so we are able to calculate the whole spectrum of shapes for different reduced volumes as shown in Figure 2. In Monte Carlo simulations, we are changing the parameter Δp. Once we obtain the equilibrium state for a given Δp, we can also calculate the average reduced volume v for that state in order to compare the thermal equilibrium and ground state calculations. As can be seen in Figure 3, the thermal equilibrium states obtained by Monte Carlo simulations are in good agreement with the ground state shapes calculated with the minimization of the bending energy.First row: snapshots of equilibrium configurations obtained by Monte Carlo simulations as described in Section2.2. Second row: ground state shapes calculated with the minimization process as described in Section 2. Pressure differences are in units kBT/lmin⁡3. All the shapes were calculated for C0=0. Parameters for first row pictures are as follows: (a) Δp=-0.01 (with v=0.949±0.008), (b) Δp=-0.05 (with v=0.769±0.003), (c) Δp=-0.073 (with v=0.66±0.01), and (d) Δp=-0.1 (with v=0.172±0.006). Shapes in the second row were calculated for the following reduced volumes: (a) v=0.95, (b) v=0.77, (c) v=0.66, and (d) v=0.17.
(a)
(b)
(c)
(d)We can see in Figure2 that the ground state oblate shapes exist only in a small region of the values for v. The fluctuations which are taken into account in Monte Carlo simulations alter the phase diagram presented in Figure 2. The regions of stability of different shapes are changed. For the values of parameters investigated in the Monte Carlo simulations, we did not observe stable oblate shapes. We cannot exclude the possibility that the stable oblate shapes can be obtained in Monte Carlo simulations, but we can speculate that the region of stability of oblate vesicles is small. The Monte Carlo evolution of the vesicle from an initial quasispherical state towards an equilibrium stomatocyte state is shown in Figure 4. As can be seen, the vesicle spends a relatively long “time” in a metastable oblate discocyte state before it reaches the equilibrium state.Figure 4
Monte Carlo evolution of the vesicle with pressure differenceΔp=-0.071 from an initial quasispherical state towards an equilibrium stomatocyte state. The reduced volume of the vesicle is shown as a function of the Monte Carlo time (measured in mcs). The average reduced volume in the equilibrium stomatocyte state is v=0.185±0.007, while in the metastable oblate discocyte state it is v=0.50±0.02. The snapshots of the initial quasispherical, discocyte, and stomatocyte states are shown for the appropriate values of the reduced volume.Near the boundary between prolate and stomatocyte equilibrium states, calculated with the Monte Carlo method, the fluctuations of the vesicle increase (as can be seen in the increased standard deviation of the reduced volume). Due to thermal fluctuations (i.e., the stochastic nature of the Monte Carlo method), the boundary between the prolate and stomatocyte vesicles cannot be determined with a high accuracy. For example, in Figure3, the equilibrium state for Δp=-0.073 is a prolate shape, while we can obtain stomatocyte equilibrium states already for Δp=-0.071 (Figure 4).We have investigated the equilibrium configuration of the nematic liquid crystal on a prolate vesicle with the reduced volumev=0.70, which was calculated in the spontaneous curvature model. The nematic ordering on a vesicle is shown in Figure 5. The topological defects are the points where λ=0 (dark red color). We have observed four topological defects, each with the charge 1/2, which means that the theorem of Poincaré [37] is satisfied. We can see that the topological defects occur in the regions of the shell with the highest curvature. Their position is strongly curvature driven, as previously described in [11, 14]. The regions of the vesicle with lower curvature have a higher degree of orientational order; therefore it is not energetically favourable for the defects to occur in those places. We can find many complex configurations of topological defects on different vesicles.Figure 5
Vector field of molecules and contour plot ofλ/λc ratio, where the condensation value λc stands for the highest possible value of λ at given t<0 reduced temperature. The calculations were performed for R/ξ0=40 and t=-0.03. The vesicle shape is prolate with the reduced volume v=0.70, calculated in the spontaneous curvature model for C0=0. A schematic representation of the vector field and the contour plot on a vesicle is on the right side.


## 4. Conclusions

We have studied the vesicles of spherical topology within the framework of the spontaneous curvature model. The shapes of the vesicles were calculated both by Monte Carlo simulations and by the minimization of the curvature energy functional. The vesicles shapes calculated by both methods are in very good agreement. In Monte Carlo simulations, the calculations were performed for a fixed value of the pressure difference. Therefore, we were not able to obtain the vesicles of the given reduced volume,v, with sufficient precision. We were able to obtain only metastable oblate shapes. We observed in the simulations that a vesicle spent a relatively long time in a metastable oblate state before it reached a stable stomatocyte state. The minimization of the curvature energy indicates that the oblate shapes are stable in a very narrow range of the reduced volume. It may be possible that the region of stability of the oblate shapes has changed due to the fluctuations or even that the oblate shapes have become only metastable.The nematic ordering on a prolate vesicle was also studied. The shape of the vesicle was calculated in the spontaneous curvature model. The net topological charge on the surfaces with the topology of a sphere equals2. The equilibrium configurations with four topological defects, each with charge 1/2, were also calculated. We were able to control the positions of the topological defects on a vesicle by changing the curvature of the vesicle. The tendency of the topological defects to accumulate in the regions of high curvature may be important for fission processes [40]. The vesicles with the nonzero spontaneous curvature have very often the shapes composed of beads connected by small necks. The accumulation of the topological defects in those necks may lead to breaking the necks and dividing a large vesicle into a few smaller vesicles.In the future research, one could calculate the nematic ordering on oblate and stomatocyte vesicles. It would be possible to use the Monte Carlo simulations in order to calculate the vesicle shapes for different value of the spontaneous curvatureC0 and compare them to the shapes calculated with the minimization procedure. The nematic ordering could also be studied on those vesicles. It would be instructive to see whether on any of those vesicles pairs of defects with opposite topological charges could be generated.



---

*Source: 373674-2014-11-30.xml*