Non-spherical shapes of capsules within a fourth-order curvature model

We minimize a discrete version of the fourth-order curvature based Landau free energy by extending Brakke's Surface Evolver. This model predicts spherical as well as non-spherical shapes with dimples, bumps and ridges to be the energy minimizers. Our results suggest that the buckling and faceting transitions, usually associated with crystalline matter, can also be an intrinsic property of non-crystalline membranes.


Introduction
The study of organized structures, like self-assembled membranes and vesicles, is an important part of soft matter physics that is relevant for chemistry and biology and for practical applications [1,2]. These systems can usually be described as two-dimensional surfaces, since their thickness is much smaller than the typical size. Thus, at the phenomenological level, the search for equilibrium shapes of such micro-and nano-structures is a problem of differential geometry of surfaces. This view is certainly appropriate for liquid membranes where the free energy is considered as a functional of the local curvatures and does not depend on the in plane deformation tensor like for crystalline membranes [3].
A simple theory for elastic shells was proposed by Sophie Germain around 1810 [4,5] in which the energy is given by the form where H is the mean curvature, K is the Gaussian curvature and the integral is over the surface S. In the case of soap films, the first term, proportional to the surface tension σ, is the dominant contribution to the free energy. For fluid membranes, instead, σ ≃ 0 because molecules can easily flow and adjust the total area of the membrane to the one corresponding to the best packing [3]. The second term in the integrand in Eq. (1) describes out-of-plane bending of elastic shells. According to the Gauss-Bonnet theorem [6] for compact surfaces, the last term dS K = 2πχ(S) is a topological invariant, with χ(S) called the Euler-Poincaré characteristic of the surface. For all surfaces topologically equivalent to spheres χ = 2 and dS K = 4π. The squared mean curvature integral W = dS H 2 is known in mathematics as Willmore functional [7], in the theory of membranes as Helfrich free energy [8], and in string theory as Polyakov action [3]. The sphere turns out to give the absolute minimum of this functional, W = 4π, among all compact surfaces [7]. This naturally explains why spherical shapes are common in the world of fluid membranes. More complicated equilibrium shapes than the spherical one, e.g. prolates, discocytes and stomatocytes [9,10], can be obtained within the Willmore functional by adding the condition of constant volume. Complicated equilibrium shapes are also observed in systems, like red blood cells or viral capsids, where other degrees of freedom rather than bending become important [3] so that the functional W and thus Eq. (1) is not sufficient. These situations are usually described in terms of crystalline membranes where these additional degrees of freedom are described by the in-plane deformation tensor.
In Ref. [11], for symmetric bolaamphiphilic fluid membranes, it was found that the interactions of the hydrophilic tails of the bolaamphiphiles with molecules of the solvent as well as entropic terms may lead to a negative coefficient A in Eq. (1). Then, symmetry allowed higher order terms should be added to stabilize the free energy that otherwise would not be bounded from below leading to the free energy functional where the minus sign is written explicitly so that from now on A is positive. The interplay between these higher order terms and the Willmore functional with negative sign leads to spontaneous bending. This model was successfully applied to describe the experimental data on deformations of spherical bolaamphiphilic vesicles in high magnetic fields [12]. Here we study the equilibrium shapes of liquid membranes with spontaneous bending based on minimization of the functional Eq. (2).
Without referring to any specific system, we minimize numerically F 4 under the constraint of constant surface (S = 4πR 2 ) and preserving the topology of the sphere (χ = 2). To this purpose we have supplemented the opensource software "Surface Evolver" [13] with new subroutines for the calculation of fourth order terms. Increasing A, for a given set of c i , we find a transition from a spherical shape to more complicated, dimpled, shapes when AR 2 /c 1 ≃ 2. This continuous transition is followed by a discontinuous one towards shapes of icosahedral symmetry with ridges and facets, when AR 2 /c 1 > 8. The presence of the coupling term between the mean curvature and the Gaussian curvature (c 2 = 0) results in intermediate shapes with bumps. We characterize these new equilibrium shapes in terms of order parameters (rotational invariants), which were introduced to describe virus capsids [14] and bond-orientational order in liquids and glasses [15]. Moreover, by comparing rotational invariants, we draw an analogy between the buckling transition found in virus capsids [14] and the transition from spherical shapes towards shapes with dimples and ridges explored in this paper.

Phenomenological models
The spontaneous-curvature model introduced by Helfrich in 1973 [8] accounts for a possible asymmetry of membranes, such as a difference in the number of molecules in each layer of bilayer vesicles. He suggested the following free energy of fluid membranes where H 0 is called the spontaneous curvature, κ is known as bending rigidity andκ is the Gaussian rigidity that affects only transitions implying a change of topology. The Helfrich model was successful in describing different phenomena, like the budding transition, discocytestomatocyte transition, et cetera [9,10]. However, for some systems it turns out to be too simple and thus insufficient to explain new experimental data. In fact the Helfrich approach accounts only for the basic degrees of freedom, such as bending, neglecting the possibility of a tilt of the molecules within a layer [16], stretching/compression of layers [17] and interaction with the environment [11]. The latter situation may result in a free energy with spontaneous bending (negative coefficient in front of H 2 ) given by Eq. (2). In general, one can imagine other mechanisms favouring spontaneous bending, like geometric constraints of packing, complex van der Waals or electrostatic intermolecular interactions.
In this paper we consider only symmetric fluid membranes, where the Helfrich free energy F H with H 0 ≡ 0 coincides with the Willmore functional W = dS H 2 . It is worth noting that F H does not depend on the size of the structure. The presence of higher order terms in F 4 (Eq. (2)), on the contrary, result in a characteristic length scale, L ∝ c i /A. To guarantee that F 4 has a minimum for real values of H and K we require the form Thus a minimum exists if To relate the coefficients of Eq. (2) to the bending rigidity κ, entering the Helfrich model Eq.
(3), we calculate In the case c 2 = 0, ∀c 1 , c 3 > 0 this expression simplifies to κ = A. In principle the parameters c i can be determined by comparing equilibrium shapes with experimental values for the deformations due to high magnetic fields as it has been done for sexithiophene vesicles [12]. However, the coefficients c i are not separately accessible. Here, we will consider them as formal parameters, satisfying Eq. (4) and study the possible equilibrium shapes and their transformations, without referring to specific systems. To find equilibrium shapes, one usually needs to solve the Euler-Lagrange equation. For the Willmore functional W only six analytic solutions are known: planes, cylinders, spheres, tori, cones, dupin cyclides [18]. Since we do not consider here topological transformations, the sphere gives the absolute minimum W = 4π for χ = 2. The equilibrium shape equation for the functional F 4 can be written following Ref. [19], which results in Here ∇ 2 stands for the Laplace-Beltrami operator, [6]. By substituting H 2 = K = 1/R 2 , where R is the radius of a sphere, into Eq. (6), one finds a sphere as solution if and only if c 1 + c 2 + c 3 = 0, which contradicts the condition (4). Therefore a sphere is not an equilibrium shape of the free energy (2). Finding analytical solutions for a highly nonlinear partial differential equation (Eq. (6)) is not likely. Alternatively, one can study equilibrium shapes by means of numerical methods. Here we adapt the "Surface Evolver" software [13,20] to this purpose.

Computational details
In the interactive program "Surface Evolver", a surface is modeled by a set of triangles, which is finite for compact surfaces. Given a triangulation of the surface we evolve it towards the shape that minimizes the total free energy. For the evolution "Surface Evolver" uses the steepest descent method, which means that at each iteration all vertices are moved along the gradients of the free energy. Then we can refine the surface by dividing each triangle into four, and repeat the procedure until the approximated surface becomes smooth. An example of the evolution, starting from the icosahedron, towards a sphere is shown in Fig. 1. Since the Euler characteristic χ does not depend on the triangulation of the surface, but only on the topology, it always holds χ = F − E + V = 2, which relates the total number of triangles (faces F ), edges (E) and vertices (V ).
The first application of "Surface Evolver" was to study the shape minimizers of the Willmore functional W , starting from polyhedra with different χ Here we aim at studying the energy minimizing shapes of the free energy (2), starting, for reasons that we explain later, from an icosahedron (see Fig.1a). Surface Evolver v.2.30 evaluates the energy terms, dS H 2 and dS K 2 . We have written two new subroutines to calculate the other two terms entering Eq. (2), namely dS H 4 and dS H 2 K [22]. For discrete surfaces at each vertex ν, the mean curvature H ν and the Gaussian curvature K ν are defined as [23,24] where S ν is a Voronoi area around ν, i.e. the area of the cell consisting of all points closer to ν than to any other vertex, ∇S ν is the gradient of S ν at ν, and i φ i,ν is the sum of all facet angles at the vertex ν of icosahedron. The definition of H ν comes from the fact that the mean (extrinsic) curvature measures the variation of the area element, displaced along the normal, divided by the corresponding change of volume. The definition of K ν comes from the Gauss-Bonnet theorem for a Voronoi region. Then, we can approximate the integrals by assigning their energy contributions to the vertices only. The integrals are calculated as a sum over all vertices ν of the curvature times the area around a vertex, which gives: Note that, for the calculation of the energy contributions, the area associated with a vertex is taken to be 1/3 of S ν , rather than S ν . This approximation simplifies the calculations for "Surface Evolver" and works best for triangles close to equilateral. The choice of icosahedron as a starting shape guarantees that, at every iteration, we are close to equilateral triangulation so that the approximation of the integrals as in Eqs. (8) and (9) holds. Moreover, among all regular polyhedra, the icosahedron has the ratio of the surface area over the enclosed volume closest to that of a sphere. It is in general convenient to start the minimization from a shape close to a sphere because the steepest descent method implemented in "Surface Evolver" finds only one stable local minimum, which is not necessarily a global one. Moreover, in nature many shapes, like viruses, are found to have an icosahedral symmetry and we will compare the predictions of our phenomenological model with the models developed for viral capsids, in particular the one discussed in [14].

Numerical results
Assuming a constant area of a surface S = 4π (R = 1), we minimize the free energy given by Eq. (2) with"Surface Evolver". In Table 1 we illustrate the equilibrium shapes as a function of A for three typical cases: i) c 1 = 1, c 2 = c 3 = 0, ii) c 1 = 1, c 2 = 0, and c 3 = 0.5, iii) c 1 = 1, c 2 = −1, c 3 = 0.5. Irrespective of the particular choice of the constants c i , spherical shapes (shapes very close to a sphere [25]) are found for negative values of A (−A dS H 2 is positive). For A ≃ 4 we see a transition towards dimpled shapes (second column in Table 1) Table 1. Summary of the equilibrium shapes for different values of the parameters in F 4 = dS{−AH 2 + c 1 H 4 + c 2 H 2 K + c 3 K 2 }. In the rows from left to right, shapes evolve from spherical towards non-spherical ones with increase of the negative contribution −A dS H 2 . The number of the column, given in brackets, corresponds to the label of the points in Fig. 3 and Fig. 4. All the shapes have a constant area and triangulation with a number of vertices V = 642.
followed by a transition towards shapes with ridges (last column in Table 1) connecting the vertices of icosahedron.
In presence of a non-zero coupling term c 2 between the mean curvature H and the Gaussian curvature K, intermediate shapes with bumps (see bottom row in Table 1) occur, favouring a positive K for c 2 < 0. Note that, the bending rigidity κ increases with A according to Eq. (5). Thus the stiffness of shapes increases in the rows, for the first and the second rows κ = A (c 2 = 0), for the third row κ = 2A.
A canonical way to characterize the shapes presented in Table 1 is to expand their radial density R(θ, φ) in terms of spherical harmonics Y lm (θ, φ), as follows, with the coefficients Q lm of the above expansion defined as In the case of triangulated surfaces we deal with a discrete radial density R(θ, φ), defined at the vertices of triangles. The coefficients Q lm were computed in Matlab using Gauss quadratures. Then, according to [15], we construct second and third order rotational invariants, as and where l l l m 1 m 2 m 3 are the Wigner 3j symbols. For a sphere the only non-zero coefficient is Q 00 = √ 4π, giving Q 0 = 4π. The shapes with icosahedral symmetry are distinguished by non-zero invariants Q l and W l for l = 0, 6, 10, 12, . . . [15]. The vanishing of invariants for other values of l was also recovered in the present calculations. For the icosahedron, only for very high refinements, namely for triangulation with number of vertices V = 10242, we found for W 6 and W 10 (notice W 0 ≡ 1) the same values as in Ref. [15]. Therefore, in the following, all the integrals in Eq. (11) are computed with V = 10242, rather than the one presented in Table 1.
We plot first Q 0 /4π, characterizing the aspherity of the shape, as a function of the adimensional parameter AR 2 /c 1 in Fig. 2. We define spherical shapes [25] when Q 0 /4π ≃ 1, implying that AR 2 /c 1 2. Non-spherical shapes correspond to AR 2 /c 1 > 2 for all three curves with different values of c 2 and c 3 . The biggest deviation from a sphere, Q 0 /4π ≈ 0.86, happens at A = −12, c 2 = 0 and c 3 = 0.5. The corresponding dimpled shape can be found in Table 1. Notice, that the volume of the equilibrium shapes decreases with increasing A in the same way as Q 0 .   Table 1. For comparison, we show the data from [14], used to describe the buckling (points b-d) and faceting (points d-f) transitions of viral capsids. Figure 3 shows the second order invariants Q l /Q 0 for l = 6, 10 plotted against each other. Starting from a spherical shape labeled by point (1) we follow the curves with the points separated by δA = 4, as in Table 1. The dimpled shapes correspond to the growth of Q 6 /Q 0 (point 2 for all curves), whereas the appearance of ridges is characterized by the increase of Q 10 /Q 0 along the curves. By adding the points from Fig. 8 in Ref. [14], we compare our equilibrium shapes for fluid membranes, with shapes of crystalline viral capsids studied within the continuum elastic theory. The governing parameter of that model, relating stretching and bending of elastic shells, is the dimension-  Fig. 4. The third order invariants W 6 and W 10 distinguish the shapes with bumps and ridges (W 6 < 0) from the dimpled shapes (W 6 > 0). The points 1-6 are labeled as for Fig. 3 and Table 1. The spherical shapes correspond to point 1 at (0, 0). The icosahedron is marked by an asterix.
less Föppl-von Kármán (FvK) number γ = Y R 2 /κ, where Y is the Young modulus. The points b-d describe the buckling transition of viral capsids, and the points d-f are associated with sharpening of the ridges at large γ [14]. In our case, the transition towards dimpled shapes may be analogous to the 'buckling', whereas the appearance of ridges leads to the increase of Q 10 /Q 0 contrary to the model of viral capsids. To find out more connections between the non-spherical equilibrium shapes and to distinguish shapes with bumps and dimples, we consider the combination involving the third order invariants W l , as shown in Fig. 4. We notice that shapes with dimples are characterized by W 6 Q 6 ≈ 0.02, whereas shapes with bumps have the same magnitude but the opposite sign of W 6 (curve with rhombuses). The sign of W 6 , which is the first non-zero third-order invariant, discriminates the shapes with dimples from the shapes with bumps and ridges. The sharpening of ridges is characterized by an increase in the invariants Q l and W l calculated for higher degree l = 10. According to Figs. 3, 4 the point 3, corresponding to A = 8, c 1 = 1, c 2 = −1, c 3 = 0.5 (see also Table 1) is the closest one to the icosahedron.

Discussion and conclusions
We have studied the equilibrium shapes of symmetric fluid membranes with a spherical topology. Assuming a fourthorder curvature model proposed in Ref. [11] we found a variety of shapes with dimples, bumps and ridges as well as quasi-spherical shapes. We noticed that similar shapes appear in the theory of elastic icosahedral shells, when studying the buckling and ridge-sharpening transitions [14,26]. Both these transitions depend only on the FvK number, which is the ratio of the stretching and bending contributions to the free energy. In our case, the competition arises between the negative quadratic term in Eq. (2) and higher order quartic terms. Our numerical analysis shows that the transition from spherical towards dimpled shapes depends only on the value of A or more likely on the dimensionless combination AR 2 /c 1 . The transitions between non-spherical shapes, such as dimples-ridges and dimples-bumps, are not determined only by the parameter AR 2 /c 1 . Both buckling and ridge-sharpening transitions occur within one order of AR 2 /c 1 whereas in the model of elastic shells the FvK number should change within four orders of magnitude [14,26].
Our calculations were done under the constraint of constant surface, and we found a decrease of the volume with increasing A, similar to the change of volume upon buckling transition [26]. It might be interesting to study the F 4 -minimizing shapes under the constraint of constant volume. With the latter constraint, non-spherical shapes appear also with the Willmore functional but those shapes are essentially different from the ones we find within our model, namely they present large deviations from spheres but no bumps, ridges or dimples.