Numerical Study of Membrane Configurations

1 Laboratory of Biophysics, Faculty of Electrical Engineering, University of Ljubljana, Tržaška 25, SI-1000 Ljubljana, Slovenia 2 Laboratory of Clinical Biophysics, Faculty of Health Sciences, University of Ljubljana, Zdravstvena 5, SI-1000 Ljubljana, Slovenia 3 Department of Physics, Faculty of Natural Sciences and Mathematics, University of Maribor, Koroška Cesta 160, SI-2000 Maribor, Slovenia 4 Jožef Stefan Institute, P.O. Box 3000, SI-1000 Ljubljana, Slovenia 5 Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 44/52, 01-224 Warsaw, Poland


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][12][13][14].On a spherical surface, the equilibrium configuration typically has Advances in Condensed Matter Physics 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 Section 2, 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.

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][18][19][20]] where  and  are bending constants,  1 and  2 are principal curvatures of a vesicle surface, and  0 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  tot is calculated by integrating (1) over the entire surface of a vesicle: where  is an infinitesimal element of the vesicle area .The Gauss-Bonnet theorem states that the last term on the righthand 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 the -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 -axis by an angle  = 2.The surface of the vesicle is constructed in this manner.To describe the vesicle profile on the  −  plane, we introduce the arc length  of the profile curve and an angle ().() is the angle of the tangent to the profile curve with the plane that is perpendicular to the axis of rotation  (Figure 1).If the function () is known, the vesicle profile can be calculated by the following parametric equations: where () and () are the coordinates of the vesicle profile in the  −  plane.We can expand the function () in the Fourier series [26]: where   is the profile length,  is the number of Fourier modes, and   are the Fourier amplitudes, which are calculated when the bending energy  tot is minimized.For closed shapes, we apply the following boundary conditions: (0) = 0, (  ) = , and (0) = (  ) = 0, which means that  0 =  in (4).The bending energy  tot is a function of Fourier amplitudes   and the profile length   ; therefore we need to do the numerical minimization of the function of many variables in order to calculate vesicle shapes [26][27][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.

Monte Carlo Simulations.
The fluid vesicle is represented by a set of  vertices that are linked by tethers (i.e., bonds) of variable length  so as to form a closed, randomly triangulated, self-avoiding surface [29,30].The lengths of the tethers can vary between a minimal value,  min , and a maximal value,  max .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  max / min = 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 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, , due to the pressure difference, Δ, 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 ( 0 = 0).The bending energy  tot 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 in Monte Carlo sweeps (mcs).One mcs consists of individual attempts to displace each of the  vertices by a random increment in the sphere with the radius , centered at the vertex, followed by    attempts to flip a randomly chosen bond.We denote   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   = 3 and / min = 0.15.
The investigated vesicle consists of  = 1447 vertices, which are connected with 3( − 2) = 4335 bonds to form   = 2( − 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, ) is allowed to vary between 4 and 8.For the bending stiffness of the vesicle, we use  = 10  , where   is the Boltzmann constant and  the absolute energy.In the following, we use  min as the unit of length and    as the unit of energy.

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 (e 1 , e 2 ) in which we can represent Q as [11,35] Q =  0 (e 1 ⊗ e 1 − e 2 ⊗ e 2 ) +   (e 1 ⊗ e 2 + e 2 ⊗ e 1 ) , (6) where  0 and   are scalar functions in the chosen coordinate system.The tensor Q can also be written in a diagonal form: 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 of  0 and   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] where  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  is negative, which means that the system is below critical temperature.The operator ∇ 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 q0 and q at random points on the surface, until we reach the equilibrium configuration.The tilde notation is used, because we scaled  0 and   in order to get a dimensionless expression for energy (see (8)).

Results and Discussion
The result of minimizing the bending energy  tot , given by ( 2), is presented in Figure 2. Three classes of vesicle shapes in their ground states (i.e., in states with the minimal  tot ) 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].
In Figure 3, 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 Δ.Once we obtain the equilibrium state for a given Δ, 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.
We can see in Figure 2 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.
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 Figure 3, the equilibrium state for Δ = −0.073 is a prolate shape, while we can obtain stomatocyte equilibrium states already for Δ = −0.071(Figure 4).
We have investigated the equilibrium configuration of the nematic liquid crystal on a prolate vesicle with the reduced volume V = 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.

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 equals 2. 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

Figure 2 :
Figure 2: Bending energy of closed vesicles in units  =  tot /(8) as a function of reduced volume V. Shapes were calculated within the spontaneous curvature model for  0 = 0.

Figure 4 :
Figure 4: Monte Carlo evolution of the vesicle with pressure difference Δ = −0.071from 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.