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.

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, 1114]. 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 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.

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, 1720] where and are bending constants, and are principal curvatures of a vesicle surface, and 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 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 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 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 . 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 is minimized. For closed shapes, we apply the following boundary conditions: , , and , which means that in (4). The bending energy 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 [2628]. 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 . The reduced volume is defined as the ratio of the vesicle volume to the volume of the sphere with the same surface area as a given vesicle.

2.2. 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, , and a maximal value, . 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 .

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 (). The bending energy 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 and .

The investigated vesicle consists of vertices, which are connected with bonds to form 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 , where is the Boltzmann constant and the absolute energy. In the following, we use as the unit of length and 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 . On the surface, we introduce a local orthonormal basis () in which we can represent as [11, 35] where and are scalar functions in the chosen coordinate system. The tensor can also be written in a diagonal form: where and are eigenvalues of eigenvectors and . Value is a measure for the orientational order and is bound to lie on an interval . The value corresponds to the maximal degree of orientational order, while the value corresponds to an isotropic state with no orientational order at all. The points on the vesicle exhibiting 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 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 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 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 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 and at random points on the surface, until we reach the equilibrium configuration. The tilde notation is used, because we scaled and in order to get a dimensionless expression for energy (see (8)).

3. Results and Discussion

The result of minimizing the bending energy , given by (2), is presented in Figure 2. Three classes of vesicle shapes in their ground states (i.e., in states with the minimal ) are shown along with their energies for different values of the reduced volume . 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 , 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 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 . 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 is a prolate shape, while we can obtain stomatocyte equilibrium states already for (Figure 4).

We have investigated the equilibrium configuration of the nematic liquid crystal on a prolate vesicle with the reduced volume , 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 (dark red color). We have observed four topological defects, each with the charge , 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.

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, , 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 . The equilibrium configurations with four topological defects, each with charge , 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 curvature 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Authors’ Contribution

Luka Mesarec, Miha Fošnarič, and Samo Penič contributed equally to this work.

Acknowledgments

WTG would like to acknowledge the support from NCN Grant no. 2012/05/B/ST3/03302. The authors wish to acknowledge the ARRS Grants J1-4109, J1-4136, J3-4108, and P2-0232.