Additional contributions to elastic energy of lipid membranes: Tilt-curvature coupling and curvature gradient
Konstantin V. Pinigin, Peter I. Kuzmin, Sergey A. Akimov, Timur R. Galimzyanov
11 Additional contributions to elastic energy of lipid membranes: Tilt-curvature coupling and curvature gradient
Konstantin V. Pinigin * , Peter I. Kuzmin, Sergey A. Akimov, Timur R. Galimzyanov * A.N. Frumkin Institute of Physical Chemistry and Electrochemistry, Russian Academy of Sciences, 31/4 Leninskiy prospekt, Moscow 119071, Russia; *corresponding authors: [email protected]; [email protected] ABSTRACT Lipid bilayer membranes under biologically relevant conditions are flexible thin laterally fluid films consisting of two unimolecular layers (monolayers) each about 2 nm thick. On spatial scales much larger than the bilayer thickness, the membrane elasticity is well determined by its shape. The classical Helfrich theory considers the membrane as an elastic two-dimensional (2D) film, which has no particular internal structure. However, various local membrane heterogeneities can result in a lipids tilt relative to the membrane surface normal. On the basis of the classical elasticity theory of 3D bodies, Hamm and Kozlov [Eur. Phys. J. E 3, 323 (2000)] derived the most general energy functional, taking into account the tilt and lipid monolayer curvature. Recently, Terzi and Deserno [J. Chem. Phys. 147, 084702 (2017)] showed that Hamm and Kozlov’s derivation was incomplete because the tilt-curvature coupling term had been missed. However, the energy functional derived by Terzi and Deserno appeared to be unstable, thereby being invalid for applications that require minimizations of the overall energy of deformations. Here, we derive a stable elastic energy functional, showing that the squared gradient of the curvature was missed in both of these works. This change in the energy functional arises from a more accurate consideration of the transverse shear deformation terms and their influence on the membrane stability. We also consider the influence of the prestress terms on the stability of the energy functional, and we show that it should be considered small and the effective Gaussian curvature should be neglected because of the stability requirements. We further generalize the theory, including the stretching-compressing deformation modes, and we provide the geometrical interpretation of the terms that were previously missed by Hamm and Kozlov. The physical consequences of the new terms are analyzed in the case of a membrane-mediated interaction of two amphipathic peptides located in the same monolayer. We also provide the expression for director fluctuations, comparing it with that obtained by Terzi and Deserno. I. INTRODUCTION The amphiphilic nature of lipid molecules leads to their self-assembly into bilayer lipid membranes under certain conditions in aqueous solutions. The membrane outer surfaces are hydrophilic, while the membrane interior is hydrophobic [1,2]. Bilayer lipid membranes are common for many living organisms, where they constitute the structural basis of plasma membranes, secretory vesicles, Golgi apparatus, lysosomes, etc. [3]. In biological membranes the lipid bilayer serves two main purposes: it acts as the weakly permeable barrier between the cell interior and the environment and plays the role of the matrix or platform for membrane proteins, which mediate various cell functions. The functioning of living organisms demands reshaping and topology changes of membrane structures. Intermediate states of such processes as endo- and exocytosis or various types of fusion and fission include strongly bent membranes. In some processes, the energy of lipid membrane deformations is believed to be the rate-limiting factor [4–6]. This motivates the development of methods for the analysis of membrane reshaping energetics. Most lipid membranes under biologically relevant conditions are flexible laterally fluid films consisting of two unimolecular layers (monolayers) each about 2 nm thick. On spatial scales much larger than the bilayer thickness, the membrane elastic energy is well determined by its shape. The classical Helfrich theory considers membrane as an elastic 2D film, which has no particular internal structure [7]. However, such a large scale approach is insufficient for an adequate description of various membrane processes, especially accompanied by alteration of membrane topology. Besides, the Helfrich’s approach is poorly applicable for the description of deformations induced by shallow membrane inclusions, phase boundaries, etc. To analyze these phenomena, the internal structure of membranes should be taken into account by the introduction of, for example, tilt deformation arising when the average direction of lipid tails deviates from the normal vector to the lipid monolayer surface. The tilt can appear at lipid domain boundaries [8,9], at the edge of membrane pores [10] or the boundary of membrane inclusions [11]. This deformation is exceptionally critical in processes involving alteration of membrane topology, such as fission and fusion [4,6,12,13]. In addition, lipid tilting can influence membrane fluctuations [14–17]. Various geometrical phenomenological theories adopting a two-dimensional approach and addressing the tilt degree of freedom have been proposed [7,18–26]. Many of these theories ([18,21,22,24–26]) focus on particular lamellar phases with 2D ordering of spontaneously tilted lipids, while others address the lamellar liquid crystalline phase (L α ), which is more relevant for biological membranes and considered in this paper. The drawback of such phenomenological theories is that they are unaware of (1) how various energy contributions originate from underlining microscopic physics of lipids and (2) whether these contributions in fact take place. Actually, lipid membranes are 3D objects rather than 2D ones as they are usually considered. One of the earliest works adapting a three-dimensional approach to lipid membranes was developed by M. Hamm and M. M. Kozlov [27]. Based on the three-dimensional elasticity theory, Hamm and Kozlov (hereafter HK) wrote the most general classical quadratic expression for the elastic energy of a monolayer and took into account the following primary features of lipid monolayers: (1) the volume incompressibility, (2) the in-plane fluidity, (3) the transverse isotropy in a flat configuration and (4) the presence of a transmonolayer stress profile (also called pressure profile or pre-stress). In HK’s theory, the energy of a monolayer includes an effective mean curvature, Gaussian curvature and tilt terms, the elastic moduli of which are expressed through microscopic membrane characteristics. Thus, HK theory yields a quadratic energy functional of the membrane, that upgrades the curvature-based Helfrich functional in two major aspects: i) the tilt field emerges; ii) the divergence of the tilt field contributes to the mean membrane curvature, leading to the new effective curvature field decoupled from the tilt. There are also other three-dimensional theories of lipid membranes, which address the tilt deformation mode: the reduction from three dimensions to two of the classical Frank theory of liquid crystals [28] and the opposing forces model [16,29,30]. However, both of them do not address the transverse isotropy of lipid monolayers, which may potentially reflect some of their elastic properties. In fact, as we show in this paper, the new substantial terms arise from transverse shear deformations, which directly reflect the transverse isotropy and therefore are not present in the aforementioned theories. HK model is widely used for the description of different membrane phenomena: a fusion [12,31–34], fission [5,6], poration [35–37], phase coexistence and phase boundary energy [38–42], and interaction of membrane inclusions with raft boundaries [43–45]. The model served as a rational basis for the development of the theory of elasticity of bolalipid membranes [46]. However, recently, M. M. Terzi and M. Deserno (hereafter TD) revisited the HK theory [47] and raised the concern about the validity of its derivation and the final expression for the elastic energy. Initial TD’s derivation [47] was further improved in Ref. [48]. TD showed that the existing form of the elastic theory is incomplete, because the coupling term between the tilt and effective curvature had been missed in the original HK energy functional. However, TD’s functional is unstable, which follows from a divergence of fluctuation spectra [47]. Thus, it is invalid for applications, which require minimizations of the overall energy of deformations and it is unclear to what physical consequences new energy contributions might lead. Therefore, the motivation of this paper is to find out the reasons for the instability, derive the most general stable energy functional suitable for various applications, and investigate the physical consequences of the new energy contributions. Here, we show that the elastic energy functional should be further appended by one more term — the squared gradient of the effective curvature. We also consider the influence of the pre-stress terms on the stability of the energy functional and demonstrate that the effective Gaussian curvature should be neglected both because of stability requirements and the exceedance of the quadratic order of smallness. In addition, we figure out the physical consequences of the new terms, using the example of a membrane-mediated interaction of two amphipathic peptides located in the same monolayer of a lipid membrane. The comparison of the theoretical expression for director fluctuations with that of TD is also provided. In Sec. II we recall the basic concepts from HK theory, including definitions of the director and tilt and some fundamental equations. We write down the expression for the elastic energy density and discuss the influence of the pre-stress on membrane stability. In Sec. III we express the terms of the energy functional via the curvature and tilt. The revisiting of the transverse shear deformation terms is also provided. In Sec. IV we discuss the order of smallness of the new terms as well as the instability issues caused by the effective Gaussian curvature. We also discuss the physical implications of new terms and their influence on membrane fluctuation spectra. II. FUNDAMENTAL EQUATIONS A. Basic notations In this work, we follow HK’s approach [27] and utilize the classical theory of elasticity [49] in order to derive the free energy functional for the lipid monolayer. We introduce a field of unit vectors n called directors. The field is defined at the dividing surface that separates hydrophobic and hydrophilic parts of lipids. Directors are assumed to be directed into the hydrophobic part of the monolayer; they characterize the average orientation of lipid tails. The measure of the director deviation from the unit normal N to the dividing surface is described by the tilt vector nT Nn N , which is parallel to the dividing surface, as T N . The vector fields of n and T are used to parameterize the elastic energy of the lipid monolayer. In the reference configuration, the monolayer is assumed to be flat with the directors being perpendicular to the flat dividing surface. We introduce a Cartesian coordinate system xyz , the xy -plane of which coincides with the dividing surface and the z -axis is directed into the hydrophobic part of the monolayer. In the reference configuration, the dividing surface is planar; its shape is given by ( , ) ( , , 0) x y x y X , and the unit normal vector by (0, 0,1) N . In this configuration, the points inside the monolayer can be parameterized as ( , , ) ( , ) . x y z x y z X X N
Now, we consider an arbitrary deformation of the monolayer. In the new configuration, the coordinates of the points inside the monolayer can be described by function ( , , ) x y z X , which is parameterized as ( , , ) ( , ) ( , , ) ( , ) x y z x y x y z x y X X n , where function ( , ) x y X describes the shape of the deformed dividing surface and ( , , ) x y z denotes the distance between the monolayer points and the dividing surface, measured along the director ( , ) x y n . It is worth mentioning that in HK’s theory ( , , ) x y z denotes the projection of this distance to the normal to the dividing surface. However, we will use our notation, as it is convenient and does not make much of a difference. Actually, we assume that any line segment ( , , ) x y z r perpendicular to the dividing surface at the point ( , , 0) x y in the reference configuration transforms upon a deformation to a curve ( , , ) x y r , and the transformation can be written as a power series in its natural parameterization ( , , ) ( , ) ( , ) ( , , 0)! n nnn x y z x y x y x yn rX X n , where ( , , ) x y z is the length of a curve part lying between the dividing surface and point ( , , ) x y z X . Hereafter, we consider only the first two terms of the expansion, i.e. ( , ) ( , ) x y x y X n . This parameterization is also used by TD and HK [27,47]. B. Elastic energy functional Denoting the Jacobian matrix of the function X by X , we can write the expression for the Green-Lagrange strain tensor U as T U X X I , where I is a unit matrix. Physically, U represents a measure of the distance change between the points inside a body. If we consider in an arbitrary point three small intersecting line segments perpendicular to each other and parallel to the coordinate axes, then the diagonal elements of U will be a measure of the length change of these segments, while non-diagonal — a measure of the angle change between the segments [50]. If no deformation occurs, then U is a zero matrix. If we assume that all elements ij u of U are small and the monolayer is transversely isotropic in the reference configuration, then according to [27] the free energy density of the volume element is given by: ( ) 12 ( ) 2 ( )2( ) 4( ) 4 ( ), ( ) l xx yy zz zzxx yy zz xx yy zzxx yy xx yy xy xz yz F u u uu u u u u uu u u u u u u (1)where i are the elastic moduli, which generally depend on z , and l , zz — the pre-existing lateral and transverse stress inside the monolayer. This equation can be derived from the general expression for the deformation energy after applying the symmetry requirements to the Taylor series expansion up to the second order in u ij [51]. In general, for physically stable systems the energy functional should be positive semi-definite up to some constant, i.e. it should be always bounded from below. Because functional (1) is quadratic, in the physically relevant cases the pre-stress terms ( ) l xx yy u u and zz zz u do not influence its stability. Thus, we can analyze the stability of the functional with zero pre-stress terms, treating the energy as the quadratic form in the variables ij u . Actually, it is more convenient to consider xx yy u u and xx yy u u as independent variables rather than xx u and yy u separately. Apparently, and must be non-negative, because both related terms ( ) 4( ) ( ) 4 xx yy xx yy xy xx yy xy u u u u u u u u and xz yz u u are non-negative. The remaining part (terms with , , and coefficients) corresponding to xx yy u u and zz u represents a quadratic form which must be positive semi-definite. Finally, the stability conditions are:
12 21 2 345 (2) С. Small deformation assumption Following the notations from [47], we denote by i the covariant derivative operator with i = 1 or i = 2 corresponding to x or y , respectively. This operator equals to simple partial derivatives when it acts on scalars and vectors. Besides, the following notations are used: i i e X , i i e X , ij i j g e e , m mk k g e e k mk m T T
T e e with , i i e e being the surface basis vectors; , mkij g g — the metric tensor and its inverse; , k m T T — the components of the tilt vector. In addition, the following equations hold: , ki i k i k ik K K
N e e N , where ik K is the curvature tensor [52]. Following the definitions given by HK and TD, we refer to ik ik i k K K T as the effective curvature tensor. In terms of basis vectors and vector Xe z , the strain tensor obtains the following form:
21 1 2 1 321 2 2 2 321 3 2 3 3
11 12 1 e e e e eU e e e e ee e e e e . (3)The deformations of the monolayer are assumed to be small such that the components u ij of the strain tensor U satisfy the condition ij u . Consider, for example, u xx at points with z = 0. Using the relation
21 11 xx u z g e , we get that g . Similarly, the consideration of u xy and u yy leads to conclusions that
12 22 , 1 1 g g . Because zz u z , we have z , i.e. z . Writing u xz at points with z = 0, we get
11 2 xz Tu z z e n T ; therefore, T ; and, similarly, T . Hence, from the assumption of ij u , the following conditions are satisfied:
11 12 22 1 2 g g g T T . (4)These terms are considered to have the first order of smallness. From the equality ji ij T T g and conditions (4), it follows that i T and i T coincide up to the first order. Therefore, i i T T T is considered to have the second order of smallness. From the equalities xz u T and yz y u T , it follows that the derivatives i are at least of the first order of smallness. In addition, we assume that the normal vector N and the tilt vector T change slowly along the dividing surface. The characteristic lengths of their change are considered to be large in comparison with the monolayer thickness h , which implies: , 1. j ji i hK h T (5) In fact, the same scaling rules as (5) were used in HK work [27]. In general, the spatial change of any variable is assumed to be small in this sense. For example: j jm i m i h hK h K . (6)To summarize the scaling rules, we treat the order of smallness as a number of fields occurring in the term regardless of the order of derivatives. For example, T , ji hK and jm i h K are considered as the first order of smallness terms, while T , ji hK and jm i h K — as the second order. The order of smallness of the squared effective curvature gradient is discussed further in detail in Sec. IV. In this theory of elasticity, we write the energy functional, accounting for terms up to the second order of smallness. D. Pre-stress terms The pre-stress terms worth more detailed discussion. It is important to note that we are aimed at the deriving of the quadratic and stable energy functional, i.e. the monolayer energy should be bounded from below. The pre-stress terms require accurate consideration to avoid instabilities: the pre-stress itself should be small in some sense. As the pre-stress has the physical dimension (energy/length ), it cannot be directly compared with dimensionless components of the strain tensor. Therefore, the comparison with other elastic moduli is more appropriate. Analyzing the stability, one should try to find the minimum of the energy; if the minimum exists, the functional is stable. Since, by its definition, the pre-stress is multiplied by a linear deformation in the initial expression for the energy (1), while the remaining part of the energy functional is quadratic in deformations, the optimal deformation corresponding to the energy minimum should be proportional to the pre-stress. This leads to the conclusion that the pre-stress should be considered small as well as deformations. As long as the pre-stress term is linear on deformations, it does not affect the stability. However, the polynomial parameterization of the original deformation fields can lead to the instability, unless the pre-stress itself is considered small. This can be illustrated by the simple example of elastic spring subjected to an external force. The energy of the spring of stiffness k which is subjected to the external force f kL can be written as:
22 20 0 ,2 2 kLkW x f x (7)where L is the undeformed length of the spring; x is its displacement (change of the spring length); xL and fkL are dimensionless deformations with being the optimal deformation that minimizes the energy;
20 0 kL is the pre-stress. The only requirement for the stability of this system is k > 0. However, writing a quadratic parameterization of the deformation, for instance, in the form
21 2 a a and keeping only the second order in we obtain the following expression:
W a a akL (8)The coefficient next to lacks positive definiteness, and there is the parameter space, which violates the energy stability, although originally it was stable at any deformations and pre-stress values. The reason for the stability loss is that at the high pre-stress the spring’s equilibrium deformation is far from being small. Although one usually assumes x to be not large, the stability condition demands the energy to be bounded from below at any deformation, no matter how large is it. Meanwhile, the neglecting of the higher than the second order in in (8) was made under the assumption of the smallness of . It is necessary to keep higher orders in for the energy functional to be stable at any values of parameters (considering k > 0). However, including higher orders makes the model nonlinear. Therefore, the equilibrium deformation should be considered as the parameter of the first order of smallness ( ); the term
22 0 a in (8) thus has the third order of smallness, and, consequently, should be neglected. In this case the expression for the spring energy yields W a akL . This expression is both quadratic and stable at any values of parameters, even for an arbitrary large pre-stress and deformation . Although the neglect of the second-order terms next to the pre-stress reduces local quantitative accuracy, it preserves qualitative system properties, which we consider to be more important for the model. Moreover, keeping such second-order terms, one transforms the pre-stress to some kind of elastic modulus, which looks unphysical. Thus, we consider only the first order in deformation terms next to l and zz . 0 III. MONOLAYER AND BILAYER ENERGY It this section we derive the lipid monolayer and bilayer energy using the basic notations and assumptions given in Sec. II. A. Incompressibility of the monolayer Lipid bilayers have a very large volume compressibility modulus, which is approximately the same as that of water [53,54]. Furthermore, a recent study shows that lipid monolayers possess a local volume incompressibility throughout their thickness [55]. Therefore, as in Refs. [27,47], we use the assumption of the local volume incompressibility of the monolayer. Mathematically, local incompressibility means that det 1 X or that yx z X X X . This condition allows expressing ( , , ) x y z via z and the deformation fields. We first consider the case when the dividing surface is locally non-stretchable. Formally, this can be expressed in the following form: det( ) 1 ij g g . In Appendix A, we give the expressions for x X e , y X e and for the corresponding vector product (see (A9)). Substituting (A9) to the incompressibility condition and solving for ( ) z , we get: G z zz z z K K K T (9)where
1, 2 i ij k mi G km i j
K K K K K ( ijij g is the Levi-Civita tensor) are the trace and determinant of the effective curvature tensor ij ij i j K K T . This equation is similar to the analogous equation derived by TD in Ref. [47], except for the coefficient standing at the term of the first power of z , which was equal to 1 in Ref. [47]. B. Energy terms There are four combinations of the strain tensor components in energy functional (1): xx yy u u , zz u , ( ) 4( ) xx yy xx yy xy u u u u u , and xz yz u u . In the following, we express them via curvature and tilt fields. The result of the consideration of the first two terms up to the first order conceptually does 1 not differ from that of HK and TD. The thorough consideration of the transverse shear combination xz yz u u yields the tilt term and two additional contributions: tilt-curvature coupling term and squared effective curvature gradient term. In fact, with the help of (3), we can write that up to the required order:
4( ) ( ) xz yz u u T . (10)This relation differs from those obtained for ( ) xz yz u u both by HK in [27] and by TD [47]. HK omitted the term . TD in Ref. [56] referred to Ref. [57] in order to explain why xz yz u u equals T . However, Reddy in Ref. [58] wrote that one of the assumptions was that “the transverse normals do not experience elongation (i.e., they are inextensible)”, which in our notations implies ( ) z z . However, from (9) it follows that the function ( ) z is more complex, and hence one may not omit the gradient of on the right-hand side of (10). This fact has already been indicated in a revised version of TD’s theory [48]. The right-hand side of Eq. (10) has a simple and illustrative geometrical interpretation as a square of a local tilt inside the monolayer. Recall that tilt field T is defined only on the dividing surface: nT Nn N , where n is the director field, and N is the unit normal to the dividing surface. However, it is possible to extend this definition to the bulk of the monolayer, introducing local tilt ( ) z T . To accomplish this, we consider planes which are parallel to the dividing surface in the reference configuration, when the monolayer is flat. Each of these planes deforms to surface ( , , ) x y z X where z is a fixed distance from the dividing surface to the chosen plane in the reference configuration. Now, if we denote the unit normal to surface ( , , ) x y z X by ( ) z N , we can introduce local tilt ( ) z T : ( ) ( )( ) z zz nT Nn N . It follows that up to the quadratic order that (see Appendix B):
4( ) ( ) xz yz u u z T . (11)In other words, the tilt field can be different at various z . For example, (0) T might be equal to zero, while ( 0) 0 z T . Such a situation is demonstrated in Fig. 1(d), where the local tilt arising at z = h (where h is the hydrophobic thickness of the monolayer) is depicted. 2 Now, we consider the term ( ) 4( ) ( ) xx yy xz yy xy u u u u u . The modulus corresponds to the lateral shear [51]. Because of the local lateral fluidity of the monolayer, we assume that , which leads to a twist modulus being equal to zero (see Appendix C). This assumption is conventional for treating fluids (for the review see [59]). The weaker assumption of the global fluidity, ( ) 0 z dz , could be made. However, the stability condition , obtained in (2), implies the equivalence of both local and global fluidity assumptions. In addition, we stress that the free energy (1) refers to an equilibrium state: the lipid director and tilt should be considered as time-averaged quantities, whereas the moduli in this free energy reflect the energy cost of a variation of these time-averaged quantities. The fluidity conditions and equality should also be interpreted in this context, albeit transverse motions of individual lipids might be hindered at small time scales. Note that HK made slightly different assumptions about the monolayer fluidity. Writing the expression for xx yy u u in terms of the combination ( ) 4( ) ( ) xx yy xz yy xy u u u u u u and the lateral strain ( , , ) x y z e e , HK obtained that the modulus corresponding to u equals l . Then they considered both cases corresponding to two fluidity assumptions: strong l and weak l dz . However, if u is considered as an independent deformation mode, one should require the corresponding modulus to be non-negative. Therefore, strong and weak assumptions are again equivalent. A more detailed discussion of the assumptions made by HK is provided in Appendix С. C. Final expression for bilayer and monolayer energy Substituting relations (3) and (10) into Eq. (1), we arrive at the following expression for the energy density of a single lipid monolayer: TT F z zK E z z K z zz zK E z z K z z K TT (12)where: ( ) ( ) ( ) l zz z z z , (13a)3 ( ) 4 ( ) ( ) 4 ( ) E z z z z , (13b) ( ) 2 ( ) T z z . (13c)Integrating the energy density with respect to z over the monolayer thickness from the lower to upper boundary of the monolayer, we get the expression for the surface energy density of a single lipid monolayer: monod m m m mgrt c e k K K k Kkk k K K T T (14)where: ( ) m k dz E z z , (15a)
0, 0 ( ) m m k K dz z z , (15b) ( ) t T k dz z , (15c) c T k dz z z , (15d) gr T k dz z z . (15e)In order to obtain the total elastic energy of the monolayer, the integration of expression (14) over the plane z = 0 can be replaced by the integration over the deformed dividing surface because of the assumption of the dividing surface local inextensibility, which implies g . Therefore, the energy density (14) is invariant under reparametrizations of the dividing surface. In comparison with the results obtained by TD [47,48], the term ( )2 gr k K in expression (14) is new. Noting that the director divergence i i n e n up to the first order equals to K , we can rewrite (14) as: monod m m m mgrt c e k K k Kkk k nT T n n (16)where ( ) n n is the gradient of the director divergence. To derive Eq. (16) we assumed that lipid heads behave like lipid tails during tilting and therefore the direction of their director coincides with that of lipid tails. It is, however, possible that they have their own director independent of the direction of lipid tails. Then, additional terms 4 corresponding to lipid heads director should be taken into account in Eq. (16). If the tilt modulus of lipid heads is high, then lipid heads stay approximately normal to the dividing surface and the corresponding director divergence may be replaced with curvature, giving rise to bare bending contributions sometimes considered in the literature [60]. D. Accounting for the lateral stretching-compression Above we considered only the deformations with a locally constant area of the dividing surface. Here, we additionally address the lateral stretching-compression deformation. We attribute monolayer deformations to an arbitrary surface inside the monolayer and derive the expressions for the elastic moduli corresponding to this new surface. The chosen surface is assumed to be parallel to the dividing surface in the reference configuration when the monolayer is flat. A Cartesian coordinate system is now chosen is such a way that the xy -plane coincides with the new chosen surface, and we denote by a local area variation of this surface. Now, vector product (A9) should be multiplied by 1 + and therefore expression (9) for z transforms to expression (A11). The expression for xx yy u u and zz u up to the first order can be written as: . xx yy zz u u u K (17)Therefore, expressions (12), (14) and (16) for the energy density of a single monolayer become: T F z zK E z zK z z z K T (18) monod m m m mgrt m cA m A mc e k K K k Kkk k K Kk k A KB k C K T TT (19) monod m m m mgrt m cA m A mc e k K k Kkk kk k AB k C nT T n nnT n (20)where the new moduli are given by: 5 ( ), A k dz E z (21a)
0, 0 ( ),
A m k dz z (21b) ( ) , A dz E z z (21c) ( ) , T B dz z z (21d) T C dz z z (21e)and the limits of the integration depend on the position of the new reference plane. To obtain the total elastic energy of the monolayer, the surface integral of expressions (19) and (20) should be performed over the plane z = 0. However, all terms in (19) and (20) are quadratic, and therefore the total energy of the monolayer can be obtained via the integration over the deformed surface without losing any contributions up to the second order. Thus, the energy becomes invariant under reparametrizations of the chosen surface. In view of conditions (2), expression (18), and consequently the energy functionals (19) and (20), are bounded from below. Note also that the coupling modulus k c coincides up to a sign with that of the stretching gradient modulus. Such a stretching gradient term was also previously derived in Ref. [61] within the assumptions similar to those made by HK [27] and was indicated as crucial to explain some experiments [62]. The new contributions to the elastic energy of lipid membranes are illustrated in Fig. 1. FIG. 1. (Color online) The illustration of the new contributions to the elastic energy of lipid membranes. For the sake of simplicity, all deformations are presented in the one-dimensional case. All shapes are exactly calculated under the incompressibility constraint. a) The undeformed lipid 6 monolayer in the reference configuration. The hydrophobic thickness of the monolayer is chosen to be 1.5 nm; the thickness of the lipid heads region is set 0.4 nm for illustrative purposes. b) The pure tilt deformation at an angle of 20°. c) The constant stretching gradient deformation with ddx nm -1 . d) The pure bending deformation of the dividing surface (no stretching and zero tilt with respect to the dividing surface) to a sine wave 0.05sin(2 x ). e) The deformation of zero curvature and varying tilt divergence: T x = 0.05sin(2 x ), where T x is the x -component of the tilt field. Panels (c), (d) and (e) show that a local tilt, which is responsible for the new contributions, can emerge inside the lipid monolayer even if there is no tilt with respect to the dividing surface or it is small as in panel (e). Such a local tilt emerges due to thickness variations which result in transverse shear deformations inside a lipid monolayer. For illustrative purposes, the local tilt in panels (c), (d) and (e) is shown for surfaces which in the reference configuration correspond to planes parallel to the dividing surface at a distance equal to the hydrophobic thickness of the monolayer. IV. DISCUSSION In this paper, we revisited the elastic energy functionals derived by HK and TD [27,47]. We demonstrated that both functionals should be appended by the additional term proportional to the squared gradient of the director divergence. In addition, the influence of the pre-stress on the stability of the energy functional was considered. It was shown that in order to ensure the stability the effective Gaussian curvature term should be omitted. Under the assumption of the inextensibility of the dividing surface, we derived that the energy functional obtained by HK should be appended by two additional terms: ( ) ( ) K T T n and ( ) ( ) K n ; the energy functional obtained by TD — by one additional term ( ) ( ) K n . Below, the key points of the present work are discussed in detail: the order of smallness of the ( ) K term, the neglection of the effective Gaussian curvature term, the physical consequences of the new energy terms and their influence on the director fluctuations, and the use of the incompressibility constraint in the energy functional. 7 A. ( ) K term There is no term proportional to ( ) K in the energy functionals obtained by HK [27] and TD [56]. In our derivation of the monolayer energy (14) the term ( )2 gr k K comes from the integral T dz z z K . Applying the scaling rules introduced in Eq. (6), we conclude that the combination ( ) h K is of the second order of smallness. Because the integration is performed over the monolayer thickness, the term ( )2 gr k K turns out to be of the same order as the terms T dz z T and T dz z z K T corresponding to the tilt and tilt-curvature coupling, respectively. Another argument to keep the term ( ) K in the energy functional is based on the stability considerations. The energy density (12) is stable in terms of boundedness from below. Indeed, in view of stability conditions (2), we have that ( ) 0 T z and
221 2 3 1 31
1( ) ( ) ( ) ( ) 2 ( ) ( ) 0( )
E z z z z z zz . Thus, the second and third terms in the energy density (12) are non-negative. As for the first term, it does not violate boundedness of the functional, being the linear term in zK in the presence of the non-negative second-order term zK . Omitting the term proportional to K ruins the boundness of the third term, leading to the instability of the whole functional. This explains the instability of the energy functional obtained by TD as discussed in Ref. [47]. Another argument is that the term T z K in (12) originally comes from the expression xz yz u u , which is always non-negative, while neglecting K breaks this property. The term ( ) K includes spatial derivatives of the second order. They appear because of the deformation parameterization in the form of ' ( , , ) ( , ) ( , , ) ( , ) x y z x y x y z x y X X n and the incompressibility assumption, which leads to expression (9) for and for its gradient: 8 z K z n . Therefore, although there are only the first-order derivatives in initial functional (1), the derivatives of the second order appear in the final relation. Actually, there are theories of elasticity that consider the energy as an explicit function of the strain gradient [63,64] and even the strain gradients of higher orders [65,66]. In our case, after the incorporation of the strain gradients up to the n -th order in energy functional (1), the order of derivatives in the final answer will be n + 1. Nevertheless, this number is always restricted, and therefore scaling rules (5), (6) cause no formal mathematical problems. Terms T , K T and ( ) K originate from transverse shear deformation. These terms were not present in previous theories since, as was already mentioned in the introduction, these theories did not address the transverse isotropy of lipid monolayers. There are no such terms in Ref. [28] because the classical Frank theory, which was originally derived for liquid crystals from other symmetry arguments, is used there. As for the opposing forces model [16,29,30], it considers only three energy contributions: the repulsion of lipid heads, surface tension and lipid stretching. The first two terms are related to xx yy u u in Eq. (1), whereas the third term corresponds to zz u . We therefore see that the opposing forces model lacks transverse shear deformations xz u and yz u . The tilt term, which is nevertheless present in the opposing forces model, corresponds to keeping the quadratic order in the pre-stress in Eq. (1) because T is present in zz u . We, however, keep only first-order terms due to stability arguments (see also Sec. IV.C), although T in the pre-stress may be regarded as a correction to the tilt modulus [67]. B. Estimation of the new moduli k c , k gr , B and C We can roughly estimate the new moduli k c and k gr , assuming that the main contribution to the tilt modulus ( ) t T k dz z comes from the hydrophobic region of the monolayer and that ( ) T z is approximately constant in this region. If the length of hydrocarbon chains in the reference undeformed state of the monolayer is l , then ( ) tT kz l , and from Eqs. (15e), (15d), (21d) and (21e) we obtain tc k lk , tgr k lk , t k lB , t k lC . Using HK’s estimation of the tilt modulus k t ≈ 50 mN/m ≈ 12 k B T /nm ( k B is Boltzmann constant, T = 300 K) [27] and common value of the 9 hydrocarbon chain length l ≈ 1.5 nm, we get k c ≈ −5 k B T , k gr ≈ 3 k B T ∙nm , B ≈ −9 k B T /nm, C ≈ 5 k B T ∙nm. C. Neglecting the effective Gaussian curvature HK [27] and TD [47] keep the effective Gaussian curvature term ij k mG km i j K K K , which formally has the second order of smallness and originates from the pre-stress term ( ) ( ) z z , where
1( ) ( ( )) ( )2 G z z K K T is the lateral strain, the expression for which is obtained in (A10). The Gaussian curvature is widely used in the elasticity theory of lipid membranes; its elastic modulus was determined both experimentally and from molecular dynamics simulations [68–70]. However, we argue that within the framework of the simplest classical theory of elasticity presented this term is mathematically and physically doubtful, as the accounting for this term leads to an excessive accuracy. Firstly, the effective Gaussian curvature term leads to the instability of the energy functional in the case of free boundary conditions for the tilt field. Underlying energy functional (1) is always bounded from below as long as (2) holds. But one violates the stability of the functional, keeping the effective Gaussian curvature (which has the second order of smallness) in the pre-stress part. This can be illustrated by the following simple example. Consider the flat monolayer with the dividing surface occupying the square region (0, ) (0, ) a a in xy -plane with a > 0. We choose a simple parameterization using a Cartesian coordinate system and assume the projections of the tilt onto x and y axes to be the functions of only y and x , respectively, i.e. ( ) x x T T y and ( ) y y T T x . This implies that the mean curvature is zero and leads to the simplified expression for the Gaussian curvature: yxG TTK y x . Now, the total energy density takes the form: ymono xd m t x y TTe k k T Ty x , (22)where m k is the effective Gaussian curvature modulus. The full energy is the integral over the (0, ) (0, ) a a region. After integration by parts of the first term, the energy is given by:
1( ( ) (0))( ( ) (0)) ( )2 m x x y y t x y k T a T T a T k T T dxdy . (23)0 Assume now that x T and y T approach the -functions with its peak on the boundary of the chosen region. For example: , , y xx y h hT e T e (24)where the multiplier in the form of the monolayer thickness h was introduced in order to keep the tilt dimensionless. We assume that projections x T and y T are positive. While ( ) x y T T dxdy h a , ( ( ) (0))( ( ) (0)) x x y y T a T T a T as . On the other hand, if we choose the tilt to approach the -function on the other side of the region, for instance, yx hT e and x ay hT e , then ( ( ) (0))( ( ) (0)) x x y y T a T T a T as . Therefore, the functional is unstable, because it lacks the lower bound regardless of the sign of m k in the case of the free boundary conditions for the tilt field. There are no such stability problems in the classical Helfrich functional [7]: G w k K c kK with k and k being the bending and Gaussian curvature moduli, respectively, because the mean and Gaussian curvatures can be written via the sum and product of the principal curvatures c , c , and one needs only to require the positive semidefiniteness of a quadratic form
21 2 1 2 k c c kc c ⇔ k k k for the stability. However, such decomposition cannot be made when the effective mean and Gaussian curvatures are considered because effective curvature tensor ij K is not symmetric. Therefore, the replacement of mean K and Gaussian G K curvatures by the effective ones ( K and G K ) in the Helfrich functional leads to the loss of boundedness, even if the functional is appended by the tilt squared. However, the addition of the twist term ij ijtw tw twi j ij k k kT K T , where ij is the Levi-Civita tensor and tw k is the twist modulus, stabilizes the functional, which in this case has the form:
22 20 twG t kw k K c kK k
T T . (25)Expression (25) is written in the form suggested by Helfrich [7], i.e. without the tilt-curvature coupling and curvature gradient. The stability of this functional is convenient to analyze in a local 1 Cartesian coordinate system with its origin in a given point of the monolayer reference surface, where the functional can be written as: tw t kw k K K c k K K K K K K k T (26)Treating it as a quadratic form in ij K plus the tilt squared, amended by the linear terms, one obtains the following stability conditions: tw tw t k k k k k k k . However, because of the lateral fluidity of the monolayer and stability conditions (2) requiring , the lateral shear modulus is assumed to be equal to zero in energy functional (1), and hence there is no twist term in the final expression for the monolayer energy (14), and therefore the effective Gaussian curvature contribution cannot be stabilized by the twist term. Secondly, an argument against the Gaussian curvature term is that the corresponding modulus ( ) m k dz z z contains only the pre-stress parameter of the system, but not any elastic moduli of underlying functional (1), which looks contradictory. Both stability and modulus arguments suggest that it is excessive to keep the second order in the pre-stress terms since it leads to the emergence of the artificial elastic modulus and instability of the initially stable monolayer. Although there is no Gaussian curvature in the monolayer energy (14), it does not imply that the monolayer deformations with a zero mean curvature and non-zero Gaussian curvature cost no energy. Formally, from (14) it follows that the energy of such deformations is zero, but actually, this contribution was just neglected from the functional as formally being of a higher order of smallness. To retain both the energetic contribution of the Gaussian curvature and the stability of the monolayer, one needs to include the term proportional to G K . In this case, m k enters the expression for the spontaneous effective Gaussian curvature. However, such a contribution to the energy is of the higher order of smallness and makes the model nonlinear in terms of principal curvatures. Keeping only the first-order term in G K is analogous to keeping only the first-order term in K and neglecting the second-order one, which obviously makes the monolayer unstable. A similar argument should be applied to the opposing forces model [16], where the functional should be amended by G K coming from the lipid chain’s conformational free energy, as this free energy acts there as a stabilizing contribution. In Ref. [71], the problem concerning the calculation of the Gaussian curvature modulus from simulations via the expression ( ) m k dz z z was stated, as this formula often yields positive values [72,73], which are beyond the stability range. In this paper, we showed that in the quadratic model both positive and negative signs of m k lead to the instability, which can be avoided only by 2 higher-order corrections. Thus, the positive values of m k found in Refs. [72,73] do not represent a problem. Actually, in the same way as ( ) dz z z is associated with the spontaneous effective curvature via the expression
0, 0 ( ) , m m k K dz z z , the expression ( ) dz z z can be viewed as the spontaneous effective Gaussian curvature, and therefore can be of either sign. This can be demonstrated by considering a constant pure bending deformation up to the fourth order under conditions of zero tilt and inextensible dividing surface. In this case, the energy density takes the following form:
22 23 2 2 2 20
1( ) ( )2 2 2 dls G G z zw z zK z K K E z zK z K K , (27)with the corresponding two-dimensional contributions being: dls m m m m m G G m m G mm G m m G m c G c m G w k K K k K k K K k Kk K K k K k KK k K k K K (28)where: ( ) , m k dz E z z (29a)
22 0, 0 ( ) , m G m k K dz z z (29b) ( ) c k dz E z z . (29c)From Eqs. (28), (28a-b), one can see that the integral ( ) dz z z plays the role similar to that of ( ) dz z z , as both of them are involved in the expressions for the corresponding spontaneous curvatures. A notion of “spontaneous warp” similar to the spontaneous Gaussian curvature was introduced by T. M. Fischer [74], where the contribution of deviatoric bending in the form
21 2 / 2 c c , reflecting the difference between principal curvatures c and c is considered. This contribution comes from lateral shear stresses, which we neglect due to the lateral fluidity of lipid monolayers. In Ref. [74] it was also indicated that such stresses relax much faster than curvature stresses. The Gaussian curvature term G G m
K K (see Eq. (28)) contains fourth powers of c , c , and therefore definitely differs from deviatoric bending, and it has a different nature, coming not from lateral shear but rather from the bending deformation as a higher-order correction to mean curvature 3 bending. G m K may be important in situations when the regions with high c c are considered. For example, the presence of membrane components with nonzero G m K may lead to their segregation to the neck region during fission, in a similar way to that proposed in Refs. [75–77]. As for the fission event, it might occur as a result of the formation of pores or other topological defects [78]. D. Physical implications In order to demonstrate the consequences of the correction of the energy functional, we consider a simple and illustrative problem, in which tilt modes still contribute significantly. We calculate the distance dependence of the energy of the interaction of two amphipathic peptides, mediated by the elastic deformations of a bilayer. At large distances between two peptides partially incorporated into the same monolayer, the induced membrane deformations are independent, and the corresponding energies are additive. However, at smaller distances, these deformations overlap and lead to the effective lateral interaction of two peptides. In general, such an interaction occurs not only between amphipathic peptides but also between other inclusions incorporated into lipid membranes [79–83]. Recently, it was shown in Ref. [84] that the interaction of two amphipathic peptides can be accurately described via a one-dimensional approach within the framework of HK’s Hamiltonian. We will use the same approach and compare the results obtained in the framework of Hamiltonian (16) with that following from the HK Hamiltonian. A shallowly incorporated amphipathic peptide is modeled as a cylinder, one side surface of which is hydrophilic and the other is hydrophobic; the axis of the cylinder is assumed to lie in the plane of the membrane. The peptide induces a tilt in the adjacent monolayer and therefore a director jump, n n n , at its boundaries. Let ( ) u H x be the shape of the dividing surface of the upper monolayer: the distance between the reference plane and the dividing surface, measured along the normal to the reference plane. We allow amphipathic peptides to rotate around their longitudinal axis, which implies the following boundary condition for the dividing surface of the adjacent monolayer: ( ( / 2)) ( ( / 2)) ( ) u u x x
H x D H x D D n n , where x is the coordinate of the peptide center; x n and x n are the director projections onto the x -axis at the left and right boundaries of the peptide; D is the width (diameter) of the peptide, which we assume to be 1.3 nm, i.e. approximately the diameter of an -helix. The parameters of the bilayer are the following: bending modulus k m = 10 k B T [85], tilt modulus k t = 12 k B T /nm [27], and zero spontaneous curvatures. Hydrophobic thickness h , which reflects the mean ordering of the 4 hydrocarbon chains of lipids [86], is set to 1.5 nm [85]. We also add lateral tension terms to the energy functional:
11 1 2 u u d dH Hdx dx and
11 1 2 l l d dH Hdx dx per unit length along the y -axis, with σ = 0.0025 k B T /nm [87] in each monolayer. We use the same qualitative estimation of the director jump as in Ref. [84], namely / 2 x x x Dn n n D h ≈ −0.8. For the sake of simplicity, we neglect stretching-compression deformation mode, assuming the elastic modulus of this mode to be large in comparison with other moduli [85]. This leaves us with two unknown moduli k c and k gr , which we vary in order to investigate their influence on the interaction of the peptides. It is important to note that k c and k gr cannot take arbitrary values because they must satisfy the Cauchy-Bunyakovsky inequality gr t c k k k , which follows from Eqs. (15c–e) and stability requirements (2). Keeping k c and k gr within this permissible condition, we solve Euler-Lagrange equations for functional (16) (see Appendix D) and find the minimum energy at various distances between two amphipathic peptides under the boundary conditions described above. Solutions of Euler-Lagrange equations for the HK Hamiltonian can be found in Ref. [39]. The results are presented in Fig. 2. Firstly, we fix k c = 0 and vary k gr (Fig. 2(a)). As we see, the HK Hamiltonian predicts the global energy minima at a distance of about 4 nm between the peptides. This energy minimum remains almost unchanged when k gr is varied: increasing the values of k gr shifts the energy profile upwards without significantly changing its shape. Predictably, the energy profiles obtained from Hamiltonian (16) approach HK’s energy profile as , 0 c gr k k . Figure 2(c) shows the energy profile with theoretically estimated values of moduli ( k c ≈ −5 k B T and k gr ≈ 3 k B T ∙nm ) compared to the profile following from the HK Hamiltonian: the depth of the energy well drops from 0.27 k B T /nm in the case of HK’s Hamiltonian to 0.15 k B T /nm in the case of Hamiltonian (16), demonstrating an almost two-fold decrease. Next, fixing k gr = 12 k B T ∙nm , we vary the coupling modulus k c in order to explore how k c effects the energy profile. Recall that, according to Eq. (15d) and stability conditions (2), k c < 0. In addition, the stability condition gr t c k k k together with k t = 12 k B T /nm implies that k c ≥ −12 k B T . Varying k c within 5 this range, we obtain the curves presented in Fig. 2(d). One can see that an increase in the absolute value of k c is accompanied by a decrease in the depth of the energy well, which eventually disappears. For example, there is no local energy minimum at k c = –11.5 k B T . In particular, it turns out that at the global minimum in HK’s case the vectors of tilt and effective curvature gradient are antiparallel at all points of the upper monolayer between the two peptides. Therefore, accounting for ( ) c k K T and ( )2 gr k K increases the energy at this point due to the fact that k c < 0. FIG. 2. (Color online) Physical consequences of the new contributions, ( ) c k K T and ( )2 gr k K , to the elastic energy of lipid membranes, illustrated by an example of a membrane-mediated interaction of two amphipathic peptides located in the same monolayer of a lipid membrane. The figures show the dependence of the elastic energy of the bilayer on the distance d between the peptides. (a) Schematic representation of two amphipathic peptides of diameter D = 1.3 nm, inducing a jump in the boundary directors, x x x x n n n n = −0.8. (b) The energy profile obtained within the framework of the HK Hamiltonian (red curve) in comparison with the profiles obtained using Eq. (16) Hamiltonian with different values of k gr and fixed k c = 0. (c) The energy profile obtained within the framework of the HK Hamiltonian (red curve) in comparison with the profile obtained using Eq. (16) Hamiltonian with theoretically estimated values of k c = –5 k B T and k gr = 3 k B T ∙nm (green curve). (d) The energy profile obtained within the framework of the HK 6 Hamiltonian (red curve) in comparison with the profiles obtained using Eq. (16) Hamiltonian with different values of k c and fixed k gr = 12 k B T ∙nm . Thus, we showed that k c and k gr significantly quantitatively alter the amphipathic peptides interaction profile, while the qualitative features of the profile remain for the wide range of parameters. This indicates that the usage of the new approach is necessary for a quantitative description of the membrane-mediated peptide interaction. For example, the existence of the potential well (global energy minimum) at the distance of 4 nm between two peptides, predicted within the framework of the HK Hamiltonian, might provide an explanation for their cooperation and pore formation in membranes, under the assumption that pores are formed in a highly stressed region between two peptides next to each other [88]. However, as shown in Fig. 2, at certain values of moduli k c and k gr in corrected functional (16), this energy minimum disappears, and therefore the cooperative assembling of two peptides can be hindered. In fact, there are experiments [89] showing a fivefold difference in leakages, induced by amphipathic peptide GALA, of fluorescent probes from bilayer vesicles formed of POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-choline) and DOPC (1,2-dioleoyl-sn-glycero-3-phosphocho-line); the difference cannot be explained by the disparity in binding affinities of the peptide to the membranes. POPC and DOPC have approximately the same elastic model parameters: the hydrophobic thickness [90,91], tilt and bending moduli [92], along with close spontaneous curvatures [93]. Thus, peptides’ interaction profiles should be similar in both lipids if they are described by HK's Hamiltonian. We hypothesize that the explanation of these experiments lies in the fact that POPC and DOPC have different values of moduli k c and k gr which in turn lead to the alteration of the GALA peptides pairwise interaction profile. As DOPC has lower leakage than POPC, the potential well in DOPC membranes should be lower than in POPC ones (see Fig. 2). This results in a shorter time spent by peptides being close to each other and, hence, in a reduced rate of the cooperative formation of pores by the peptides. We point out that, in the case of nonzero tilt-curvature coupling modulus k c and zero curvature gradient k gr (the TD Hamiltonian), the variational problem leading to the Euler-Lagrange equations is ill-posed. In this case, one of the roots of the characteristic polynomial of Eq. (D2) is purely imaginary, giving rise to oscillating deformation modes, which are energetically unbounded from below and because of which the energy cannot be integrated over an infinite interval. However, it is possible to reformulate the problem, considering the peptides in a box with periodic boundary conditions. Nevertheless, as expected, if k c ≠ 0 and k gr = 0, the membrane is unstable: the oscillations are energetically unbounded from below. To be specific, we provide a particular example: staying 7 within the one-dimensional case, we consider a box of length 30 nm with periodic boundary conditions and two peptides, the distance between which is equal to 14.2 nm (the exact values of the chosen parameters are not important). The analysis shows that the energy of this configuration per unit length along the axis of the translational symmetry scales as
20 0 aA bA c , where A is the amplitude of the oscillations, given in nanometers, and a , b , c are known coefficients ( a ≈ −140 k B T /nm , b ≈ 6.7 × 10 −3 k B T /nm , c ≈ 1.4 k B T /nm at k c = −5 k B T , k gr = 0 k B T ∙nm , and the rest of the membrane elastic parameters being the same as used in Fig. 2). Since a > 0, it means that the oscillations are unbounded from below: the energy approaches minus infinity as A increases (see Fig. 3). If we now consider the case of two amphipathic peptides in an infinite membrane, separated by the same distance of 14.2 nm as in Fig. 3, the energy nevertheless does not have a lower bound. As has already been said, the arising oscillations cannot be integrated over an infinite interval. However, if we set the amplitude of the oscillations to zero at infinity as a boundary condition, the oscillations in the region between two peptides still exist and the energy scales as
20 0 a A b A c , where a −85 k B T /nm , b −3 k B T /nm , c k B T /nm, and therefore, since a , when A approaches infinity, the energy tends to minus infinity, implying that the membrane is unstable. Thus, since the energy does not have a lower bound, the variational problem is meaningless in the case of k c ≠ 0 and k gr = 0. FIG. 3. The lipid membrane is unstable if tilt-curvature coupling modulus k c is nonzero, while curvature gradient modulus k gr is zero. The instability is caused by membrane oscillations which are energetically unbounded from below. The figure shows two amphipathic peptides (black ellipses) embedded in a box of length 30 nm with periodic boundary conditions. The distance between the peptides is 14.2 nm, and the analytically calculated shapes of the dividing surfaces of the upper and lower monolayers are shown. a) k с = −5 k B T , k gr = 3 k B T ∙nm (as estimated in Sec. IV.B). 8 b) k с = −5 k B T , k gr = 0 k B T ∙nm , amplitude A of the oscillations is 0.5 nm. The rest of the membrane elastic parameters are the same as in Fig. 2. In panel (a) the membrane is stable, while in panel (b) the membrane is unstable due to the arising oscillations: energy w of the membrane inside the box scales as
20 0 w aA bA c , where a ≈ −140 k B T /nm , b ≈ 6.7 × 10 −3 k B T /nm , c ≈ 1.4 k B T /nm (the energy is considered per unit length along the axis of symmetry), implying that when oscillations amplitude A increases, the energy of the membrane inside the box approaches minus infinity. In this paper, we have concentrated only on the membrane-mediated interaction of amphipathic peptides and compared it with previous results. It however does not limit possible implications of new energy terms when applied to other problems, such as the interaction of transmembrane proteins [79–82] or hydrophobic inclusions [83]. In addition, membrane inclusions may induce membrane softening [94–100]. Among the mechanisms of such a phenomenon might be the concentration-curvature coupling [101] or the anisotropy of inclusions [100,102]. Also, the tilt deformation mode might lead to membrane softening [60]. Given that the new energy terms as well as the tilt mode correspond to the transverse shear deformation, they may also modulate this effect. E. Fluctuation spectra A fluctuation analysis is widely used in molecular dynamics simulations for estimations of membrane elastic moduli [14,15,106–112,16,17,47,48,73,103–105]. In this section, we provide theoretical results for the director fluctuations analyzed within the framework of the new elastic energy functional and compare them with those obtained by TD [47] and used in Ref. [48] to fit the fluctuation data. But before that, we would like to carefully consider the approach which is used to analyze tilt and director fluctuations. Initially, such an analysis was introduced in Ref. [14] for tilt fluctuations and in Ref. [17] for director fluctuations. The approach consists in prescribing Boltzmann probability distribution exp( ( ) / ) B P F k T g to find a system in a state g , where F is the elastic free energy of a bilayer. In Ref. [14], the HK Hamiltonian without the twist term is used for F , whereas in Ref. [17] F is equivalent to the HK Hamiltonian in the Monge gauge amended by the twist term. We recall that the HK Hamiltonian as well as Hamiltonian (20) is the equilibrium free energy while director n and tilt T are timed-averaged quantities. But in simulation analyses, instantaneous values of n and T are measured. However, it is not apparent whether the dynamics of individual lipid molecules is in compliance with the same free energy F , which reflects time-averaged equilibrium properties. In other words, it is not obvious whether instantaneous values of various quantities, such as the director or tilt, 9 measured in simulations can replace the mean ones in the energy (20) or HK Hamiltonian to predict the energy of a given state and then the spectrum. Following [47] and [16], we analyze fluctuations in terms of longitudinal and transverse components of vectors in Fourier space. Using incompressibility condition (A11) up to the first order, we get / 2 u u u H M h h h n and / 2 l l l M H h h h n , where indices u and l correspond to the upper and lower monolayers, respectively, H and M are shapes of the dividing and monolayer interface surfaces, respectively, is stretching of the dividing surfaces, h is the hydrophobic thickness of undeformed monolayers. We substitute these conditions into (20) and obtain the coupling matrix for the functions { / 2 u l H H H , _ u l H H H h , M , / 2 u l n n n , n n n / 2 u l } in Fourier space, from which we get for longitudinal components of n – a rather bulky expression, see Appendix E, Eq. (E1). To compare with the case of no stretching, considered by TD [47], we take the limit A k and get:
22 B2 2 tq gr t c m t k k Tq n k k k q k k , (30)where q is the wavenumber modulus. Note that gr t c k k k due to stability requirements, and therefore q q n is a monotonically decreasing function of q . More general expression (E1) is also a monotonically decreasing function of q due to stability requirements (2) (see the proof in Appendix E). In contrast to these predictions, the expression q q n increases in simulations at large q [17,48,73,113–115], while at small q , it behaves differently depending on algorithms of fluctuation analyses: in Refs. [17,73,113,114] it is constant, in Ref. [115] it decreases and in Ref. [48] it is variable. Simulations data for q q n spectrum were fitted by TD in Ref. [48]. However, this was achieved at the cost of the unstable energy functional. TD’s expression for longitudinal fluctuations can be obtained from (30) by setting k gr = 0, since there is no term proportional to K in their functional [47]. The same is true for the functional derived in Ref. [48]. Thus, one can see that the spectrum predicted by TD monotonically increases up to the point of divergence / m t c q k k k
0 and becomes inapplicable at higher q . This divergence is the very feature that allows fitting simulation spectra. Thus, so far the longitudinal spectrum of the director has been fitted only with the divergent expression in Ref. [48], but it seems rather artificial since it is achieved by using the unstable energy functional derived by neglecting the second-order term ( ) K . In Fig. 4, we demonstrate the comparision between the fluctuation spectrum of q q n as predicted by the HK Hamiltonian (Eq. (14) with k c = k gr = 0), TD Hamiltonian (Eq. (14) with k c ≠ 0, k gr = 0) and the Hamiltonian presented in this work (Eq. (14) with k c ≠ 0 and k gr ≠ 0); for k c and k gr we use the values estimated in Sec. IV.B ( k c ≈ −5 k B T , k gr ≈ 3 k B T ∙nm , T = 300 K), and for k m and k t we use the same values as in Sec. IV.D ( k m = 10 k B T , k t = 12 k B T /nm , T = 300 K). If k gr = 0, then at high q -values q n becomes negative which implies that the coupling matrix is no longer positive definite, and this causes the instability. The instability can be simply illustrated in the case of a one-dimensional infinite monolayer subjected to periodic boundary conditions with period L . In this case, the TD Hamiltonian, per unit length along the axis of the translational symmetry, is m t c dx k n k t k t n , where t and n are projections of the tilt and director onto the x -axis which is parallel to the undeformed monolayer. This expression can be rewritten as
22 2 2 cm t ct t kdx k n k t k n nk k . If we now choose t c k t k n and sin( ) n n qx , where q m L , the energy of the membrane over period L becomes / 4 t m t c n L k q k k k q . This energy is negative at / m t c q q k k k and approaches minus infinity as q , implying the instability of the membrane. Because the TD Hamiltonian is unbounded from below, i.e. it is unstable, solving Euler-Lagrange equations for such a Hamiltonian is meaningless. We point out that the application of the lateral tension, / 2 u H + / 2 l H , does not help to recover the stability of the membrane. In fact, Eq. (30) in the case of k gr = 0 and nonzero , becomes:
22 B2 2 4 2 tq c m t m c t q k k Tq n k q k k k k q k , which is nevertheless negative at large q : the numerator is always positive, as tilt modulus k t and tension are positive, 1 while the denominator is negative at large q because, although m t m c t k k k k q k is always positive ( k m > 0 and c k ), q to the fourth power has negative coefficient c k . An alternative explanation of such discrepancies involves a microscopic noise [17,113], which is proposed as the cause of the spectrum increase at large q in simulations. At the same time, at small q , the spectra of longitudinal components behave differently depending on algorithms of the analysis of fluctuations, which can be consistent with the predicted monotonic decrease of spectrum (E1). In order to assess the dependence of function (E1) on q , we use the values of moduli, theoretically estimated in Sec. IV.B: k c ≈ −5 k B T , k gr ≈ 3 k B T ∙nm , B ≈ −9 k B T /nm, C ≈ 5 k B T ∙nm, and take k m = 10 k B T , k A = 30 k B T /nm [85] ( T = 300 K) and k t = 12 k B T /nm [27] . At these values of the moduli, function (E1) depends only on bending-compression coupling modulus A , which lies within the interval , 17,17 m A m A k k k k k B T /nm due to stability requirements. It turns out that the decrease of the function within the range q is less than 20 % at a rather large range of the values of A (from –9 to 15 k B T /nm). Thus, the decrease of q q n is expected to be not large at small values of q , when the microscopic noise is not significant. We also point out that
22 B0 2 q q m A k Tq n k A k , which coincides with a well-known result of Ref. [17] up to the definition of the bending-stretching coupling modulus A . 2 FIG. 4. (Color online) The fluctuation spectrum of q q n as predicted by the HK ( c gr k k , red curve), TD ( c k , gr k , cyan curve) and Eq. (14) ( c k , gr k , black curve) Hamiltonians, where c k and gr k are the tilt-curvature coupling and curvature gradient modulus, respectively. The spectrum is constant in the case of the HK Hamiltonian, monotonically increasing and monotonically decreasing in the case of the TD and Eq. (14) Hamiltonians, respectively. Although the TD Hamiltonian rightly reflects an upward trend observed in simulations at high q -values, it is divergent (the corresponding asymptote is depicted as a dashed line). Fixing the divergence leads to the Eq. (14) Hamiltonian, which is monotonically decreasing. For the transverse component of the director fluctuation, we obtain a constant spectrum q t k Tn k , whereas in simulations this spectrum is a decreasing function of q [16,17,48,73,113]. To fix this deviation from the theoretical prediction, in Refs. [16,48] the twist contribution tw k T , derived by HK in Ref. [27], was added to the energy functionals, which leads to the altered expression: / 2 q t tw n k T k k q . However, in Appendix C, we showed that HK’s derivation is incomplete. More importantly, in Sec. III.B we showed that tw k should be equal to zero due to the equivalence of strong and weak lateral fluidity conditions. Therefore, the discrepancy between the theory and simulations remains. 3 To tackle this issue, we recall the underlying assumptions made to derive energy functional (20). First of all, expression (20) is the equilibrium free energy. It means that , , T n should be considered as time-averaged quantities. Secondly, the corresponding moduli in (20) reflect the energy cost of variations of these time-averaged quantities. Using these two assumptions along with the lateral fluidity of monolayers, we concluded that lateral shear modulus and, therefore, tw k are equal to zero. However, although they are thermodynamically liquid, at small time scales ~ 0.1 ns, which appear in fluctuation analyses [17], lateral modes of lipid motion might be hindered. This can impact fluctuations, inducing nontrivial dependence of q n on q . In view of the above, we suggest that the membrane fluctuations, which include the dynamics of individual lipids, cannot be described in detail within the framework of equilibrium free energy functional (20). F. Incompressibility assumption in the elastic energy functional An additional subtle point is worth being considering. The energy functional in (1) is defined on the set of all possible symmetric tensors U , whereas the incompressibility condition restricts U to the set satisfying the condition ' det( ) 1 2 1 X U 1 . Because of this, the energy also must be defined over this set, but then the Taylor series cannot be written in the form of (1), as the derivatives of F with respect to ij u are not properly defined. This issue can be overcome by extending the energy function to the set of all symmetric tensors [116]. Another approach is to express, for example, zz u through other components of U from the incompressibility condition, which implies that zz xx yy u u u up to the first order. From the lateral fluidity assumption, it follows that the energy density up to the second order depends only on the combination xx yy u u rather than on xx u and yy u separately. Therefore, the Taylor series can be written as
1( ) ( ) 4 ( )2 zz zz xz yz
F z u E z u z u u . Nevertheless, the final answer for the two-dimensional energy density of the monolayer does not change. ACKNOWLEDGMENTS 4 The work was supported by the Ministry of Science and Higher Education of the Russian Federation and by the grant of the President of the Russian Federation MK-3119.2019.4 . TRG was supported by the Russian Science Foundation grant i i e X , the vector product e e , i e and the lateral strain e e . We denote by i the covariant derivative operator with i = 1 or i = 2 corresponding to x or y , respectively. This operator equals to simple partial derivatives when it acts on scalars and vectors. Besides, the following notations are used: i i e X , i i e X , ij i j g e e , m mk k g e e , k mk m T T
T e e with , i i e e being the surface basis vectors; , mkij g g — the metric tensor and its inverse; , k m T T — the components of the tilt vector. In addition, the following equations hold: , ki i k i k ik K K
N e e N , where ik K is the curvature tensor [52]. Following the definitions given by HK and TD, we refer to ik ik i k K K T as the effective curvature tensor. The expression for i e was given in Ref. [47]: [ ] [ ] k kk k ki i i i k ik i i k i K T T K A B e e N e N , (A1)where we defined ki A and i B as: [ ],[ ]. k kk ki i i iki ik i A K TB T K (A2)In the case of g , the vector product e e can be written as: ' ' ij i j e e e e , (A3)where ij is the Levi-Civita tensor. Defined as ijij g , where ij is the Levi-Civita symbol, the Levi-Civita tensor obeys the following identities: , , , 2, . j ij i ij mn m n n mi j ij i ij kj k ij ij i j i j e e N N e e (A4)Using these properties, we can write: [ ] [ ]( ) . l m l m mi j i l i j m j i j l m i j ml l m m lp pi j l i j lm i j mp j i lpl m m m pi j lm j i i j mp A B A B A A B AA B A A B A B AA A B A B A e e e N e N e e N ee N N e eN e (A5)6 Hence, ' ' l m m mij ij ij pi j i j lm j i i j mpl m m mij i j j i pi j lm m p m p j i i jl m m m m mij pi j lm p m m p m p p ml m m mij pi j lm p m m p
A A B A B AA A B A B AA A B A B A B A B AA A B A B A e e N eN eN eN e (A6)The terms in front of N and p e in Eq. (A6) can be rewritten, up to the required order, as: l mij i j lm G A A K K T (A7) ( ) ([ ][2 ][ ][ ]) , m m k kp m m p pk p kmk m mmk m p p p p B A B A T K K TT K K T (A8)where i i e is the surface gradient operator. Therefore,
21 2 ( ) ij pi j G p
K K e e e e T N e . (A9)By definition, e e . Therefore, from (A9) we get:
1( ) ( ) .2 G K K T (A10)To obtain the expression for ( ) z in the case of nonzero stretching , we note that vector product (A9) should be multiplied by . Substituting this vector product to incompressibility condition e e e and solving for ( ) z up to the quadratic order gives the following expression: G z zz z z K z K K K T (A11) 7 APPENDIX B: EXPRESSION FOR A LOCAL TILT T ( z ) Here, we derive the expression for a local tilt T ( z ) inside the monolayer in order to provide a geometrical interpretation of Eq. (10) for the transverse shears:
4( ) ( ) xz yz u u T . We recall that the tilt vector T is defined as nT Nn N with N being the unit normal vector to the surface ' ( , , 0) x y X and n — the director. By definition, ( ) ( )( ) z zz nT Nn N , where N ( z ) is the normal vector to the surface ' ( , , ) x y z X with z being fixed. Therefore, two vectors T and T ( z ) coincide, if z = 0. Using the definition of the surface normal, we obtain: ' ' ' ' ' ' ' ' ij i j gz g g g g N e e e e e ee e (B1)Expressions for ' ' ij i j e e and ( ) are given in (A9) and (A10), respectively. Thus, up to the required order:
1( ) [ 1 ) ]1 ( )1(1 ( ) ) .2 ( pG ppp z K K N T N eN e (B2)From the definition of the director we get: pp z z zz z n T N T T N TT N e T (B3)and hence:
1( ) [( ) ( ) ]2 pp z z T N T T T e . (B4)From this equation, it follows that up to the second order:
1( ) ( ) ( )4 z z K z K T T T T . (B5) 8 APPENDIX C: HK'S FLUIDITY ASSUMPTIONS In this appendix, we discuss the fluidity assumptions made by HK, showing how the combination xx yy u u can be expressed via the combination ( ) 4( ) xx yy xz yy xy u u u u u u and the lateral strain e e . Using the definition of the strain tensor (3), we can rewrite
4( ) xx yy xy u u u as:
4( ) ( ) 1(1 ) 1 2 2( ). xx yy xy xx yy u u u u u e e e e e ee e (C1)Hence, if we denote ( ) 4( ) xx yy xz yy xy u u u u u by u , we have: ( 1) ( 1) . xx yy u u u (C2)From this equation, the combination xx yy u u can be evaluated as: xx yy u u u . (C3)From (C3) and energy equation (1), it follows that the modulus corresponding to u is l . Then HK considered separately the conditions of the local fluidity l and the global fluidity l dz . However, if u is considered as an independent deformation mode, one should require l for the stability of the functional. Hence, the condition l dz is equivalent to l . Furthermore, in HK's equation, ( ) 4 G u z K K T , (where ij i j T T with ij being the Levi-Civita tensor) the second-order terms, which generally should be included, are missed. Indeed, we can write u as: 9
1( ) 4( ) ( ) 4 4( ) 441 2414 2 . xx yy xx yy xy xx yy xy xx yy xy u u u u u u u u u u u ug g z K K g z K Kg g gz g g K K g K Kz K K K K (C4)From this equation, it follows that the terms multiplied by z to the zeroth power and to the first power also have the second order of smallness. 0 APPENDIX D: EULER-LAGRANGE EQUATIONS FOR THE ONE-DIMENSIONAL CASE In this appendix, we derive and solve Euler-Lagrange equations for functional (16) for the one-dimensional case. We consider a bilayer which is translationally symmetric along the y -axis. We choose the direction of the x -axis along the bilayer plane, and of the z -axis — perpendicular to it. Let ( ) u u H H x and ( ) l l H H x denote the shapes of the dividing surfaces of the upper and lower monolayers, ( ) M M x is the shape of the membrane mid-surface, ( ) u u n n x and ( ) l l n n x are projections of the upper and lower director fields ( u n and l n ) onto the x -axis. We now note that u u d ndx n and l l d ndx n . Incompressibility condition (9) up to the first order can be written as
1( ) ( ) 2 u u dH x M x h h ndx and
1( ) ( ) 2 l l dH x M x h h ndx , where h is the thickness of the hydrophobic part of the monolayer. We also add lateral tension terms to the energy functional:
11 1 2 u u d dH Hdx dx and
11 1 2 l l d dH Hdx dx per unit length along the y -axis. The tilt fields of the upper and lower monolayers are
22 2 u u u d dT N M h ndx dx and
22 2 l l l d dT n M h ndx dx . We write the bilayer energy per unit length along the y -axis as the sum of the respective energy of both monolayers and obtain the following Euler-Lagrange equations:
22 42 2 2 2 42 42 32 2 2 344 2 2 3242 32 22 c c grcc gr u cc hd dh l l n h l h l sh n ndx dxd d dM h l l n Mdx dx dxdh l h l sh n ndxd d dn s M h l sh ndx d s x lx dh (D1)where u b n n n , u b n n n , / m t l k k , / c c t l k k , / gr gr t l k k , / t s k . From here, we obtain the following equation for M : 1 c uc gr gr c gr gr l s l s ld d s dM M Mdx l l s l dx l l s l dx . (D2)The solution of this equation has the form: d x d x d x d x M C e C e C e C e C x C , (D3)the constant parameters d i ( i = 1, ..., 4) of which are too bulky to be presented here. Then, it is straightforward to get solutions for n and n by substituting (D3) into (D1). To apply the obtained solutions to the case of two interacting amphipathic peptides, we subdivide the membrane into five regions: the bilayer between the peptides, two monolayers below the peptides, and two bilayers lying between the left and right peptide and minus and plus infinity, respectively. These regions are connected in such a way as to provide the continuity of the director field and dividing surfaces. The obtained solutions can be used for the bilayer regions. The Euler-Lagrange equations for the individual monolayers are: l c lc l cc gr u l l d d dn s M h l sh ndx dx dxd d dM h l l n h l sh Mdx dx dxdh l h l sh n ndx (D4)Using linear transformations, we obtain the equation for M , which coincides with (D2). Substituting corresponding solution (D3) into (D4), it is straightforward to obtain the expression for l n . 2 APPENDIX E: Proof that q q n is a monotonically decreasing function of q In this section, we prove that the general expression for q q n is a monotonically decreasing function of q , where q is the wavenumber modulus. To accomplish this, we will show that the derivative of this function with respect to q is always negative, which follows from the stability requirements. The expression for q q n is:
22 2 2 2 2 3 4B 2 2 2 2 [ ] [ ], q c t A t c gr t gr c cA gr m c t c m A c t m A q nB k k q k k k T C k k k B k BCk k qAC k k k k k ABk k B k k q k k k A (E1)Firstly, we rewrite (E1) in the following way:
222 2 04 24 2 0 ( ) q N q Nq n f q D q D q D , (E2)where , , , , N N D D D denote the corresponding coefficients in (E1). Then, taking the derivative of (E2) with respect to q , we get: u D N u D N u D N D Ndf qdq D u D u D , (E3)where u q . The discriminant of the quadratic polynomial in brackets in the numerator of (E2), after substitution of N , N , D , D , D in terms of moduli, can be rewritten in the following form: Disc 2 2 t c t a c a t k AB Ak k Bk k Ck k , (E4)where is the determinant of the quadratic form in , , K T , which corresponds to the transverse shear deformation: T grt m c c z z z K dzkk k K K B k C K TT T T (E5)3 As ( ) 0 T z due to the stability requirements, this integral is always positive as well as the corresponding quadratic form in T , K and in the second line of (E5) (we do not consider the unphysical case of ( ) 0 T z ). Therefore, is positive and discriminant (E4) is negative. At the same time, D N — the coefficient of u to the second power in the numerator of (E3) — equals c t B k k , as c t B k k due to it being one of the principal minors of quadratic form (E5). Hence, ( ) 0 df qdq , and, consequently, f is a monotonically decreasing function of q . 4 REFERENCES [1] M. Montal and P. Mueller, Proc. Natl. Acad. Sci. , 3561 (1972). [2] F. Szoka and D. Papahadjopoulos, Annu. Rev. Biophys. Bioeng. , 467 (1980). [3] L. D. Nelson and M. M. Cox, Lehninger Principles of Biochemistry (W. H. Freeman, San Francisco, 2004). [4] P. I. Kuzmin, J. Zimmerberg, Y. A. Chizmadzhev, and F. S. Cohen, Proc. Natl. Acad. Sci. , 7235 (2001). [5] P. V. Bashkirov, S. A. Akimov, A. I. Evseev, S. L. Schmid, J. Zimmerberg, and V. A. Frolov, Cell , 1276 (2008). [6] A. V. Shnyrova, P. V. Bashkirov, S. A. Akimov, T. J. Pucadyil, J. Zimmerberg, S. L. Schmid, and V. A. Frolov, Science. , 1433 (2013). [7] W. Helfrich, Z. Naturforsch. C , 693 (1973). [8] S. L. Veatch and S. L. Keller, Biophys. J. , 3074 (2003). [9] C. Dietrich, L. A. Bagatolli, Z. N. Volovyk, N. L. Thompson, M. Levi, K. Jacobson, and E. Gratton, Biophys. J. , 1417 (2001). [10] I. G. Abidor, V. B. Arakelyan, L. V. Chernomordik, Y. A. Chizmadzhev, V. F. Pastushenko, and M. P. Tarasevich, J. Electroanal. Chem. , 37 (1979). [11] E. M. Pasini, M. Kirkegaard, P. Mortensen, H. U. Lutz, A. W. Thomas, and M. Mann, Blood , 791 (2006). [12] Y. Kozlovsky and M. M. Kozlov, Biophys. J. , 882 (2002). [13] V. S. Markin and J. P. Albanesi, Biophys. J. , 693 (2002). [14] E. R. May, A. Narang, and D. I. Kopelevich, Phys. Rev. E , 021913 (2007). [15] E. R. May, A. Narang, and D. I. Kopelevich, Mol. Simul. , 787 (2007). [16] M. C. Watson, E. S. Penev, P. M. Welch, and F. L. Brown, J. Chem. Phys. , 244701 (2011). [17] M. C. Watson, E. G. Brandt, P. M. Welch, and F. L. Brown, Phys. Rev. Lett. , 028102 (2012). [18] P. Nelson and T. Powers, J. Phys. II Fr. , 1535 (1993). [19] M. M. Müller, M. Deserno, and J. Guven, Phys. Rev. E , 061407 (2005). [20] J. B. Fournier, Eur. Phys. J. B , 261 (1999). [21] F. C. MacKintosh and T. C. Lubensky, Phys. Rev. Lett. , 1169 (1991). [22] U. Seifert, J. Shillcock, and P. Nelson, Phys. Rev. Lett. , 5237 (1996). [23] J. V. Selinger and J. M. Schnur, Phys. Rev. Lett. , 4091 (1993). [24] H. Jiang, G. Huber, R. A. Pelcovits, and T. R. Powers, Phys. Rev. E , 031908 (2007). [25] T. C. Lubensky and J. Prost, J. Phys. II , 371 (1992). 5 [26] W. Helfrich and J. Prost, Phys. Rev. A , 3065 (1988). [27] M. Hamm and M. M. Kozlov, Eur. Phys. J. E , 323 (2000). [28] D. J. Steigmannn, Int. J. Non. Linear. Mech. , 61 (2013). [29] S. May and A. Ben-Shaul, Biophys. J. , 751 (1999). [30] S. May, Eur. Biophys. J. , 17 (2000). [31] R. J. Molotkovsky, T. R. Galimzyanov, I. Jiménez-Munguía, K. V. Pavlov, O. V. Batishchev, and S. A. Akimov, Int. J. Mol. Sci. , 2598 (2017). [32] R. J. Molotkovsky and T. R. Galimzyanov, Biol. Membr. , 90 (2019). [33] R. J. Molotkovsky, V. V. Alexandrova, T. R. Galimzyanov, I. Jiménez-Munguía, K. V. Pavlov, O. V. Batishchev, and S. A. Akimov, Int. J. Mol. Sci. , 1483 (2018). [34] S. A. Akimov, R. J. Molotkovsky, T. R. Galimzyanov, A. V. Radaev, L. A. Shilova, P. I. Kuzmin, O. V. Batishchev, G. F. Voronina, and Y. A. Chizmadzhev, Biol. Membr. , 14 (2014). [35] S. A. Akimov, P. E. Volynsky, T. R. Galimzyanov, P. I. Kuzmin, K. V. Pavlov, and O. V. Batishchev, Sci. Rep. , 12509 (2017). [36] S. A. Akimov, P. E. Volynsky, T. R. Galimzyanov, P. I. Kuzmin, K. V. Pavlov, and O. V. Batishchev, Sci. Rep. , 12152 (2017). [37] S. A. Akimov, V. V. Aleksandrova, T. R. Galimzyanov, P. V. Bashkirov, and O. V. Batishchev, Biol. Membr. , 270 (2017). [38] T. R. Galimzyanov and S. A. Akimov, JETP Lett. , 463 (2011). [39] T. R. Galimzyanov, R. J. Molotkovsky, M. E. Bozdaganyan, F. S. Cohen, P. Pohl, and S. A. Akimov, Phys. Rev. Lett. , 088101 (2015). [40] G. Staneva, D. S. Osipenko, T. R. Galimzyanov, K. V. Pavlov, and S. A. Akimov, Langmuir , 1591 (2016). [41] J. J. Williamson and P. D. Olmsted, Phys. Rev. Lett. , 079801 (2016). [42] T. R. Galimzyanov, R. J. Molotkovsky, F. S. Cohen, P. Pohl, and S. A. Akimov, Phys. Rev. Lett. , 1 (2016). [43] T. R. Galimzyanov, A. S. Lyushnyak, V. V. Aleksandrova, L. A. Shilova, I. I. Mikhalyov, I. M. Molotkovskaya, S. A. Akimov, and O. V. Batishchev, Langmuir , 3517 (2017). [44] K. V. Pinigin, O. V. Kondrashov, I. Jiménez-Munguía, V. V. Alexandrova, O. V. Batishchev, T. R. Galimzyanov, and S. A. Akimov, Sci. Rep. , 4087 (2020). [45] K. V. Pinigin, M. V. Volovik, O. V. Batishchev, and S. A. Akimov, Biol. Membr. , 337 (2020). [46] T. R. Galimzyanov, P. I. Kuzmin, P. Pohl, and S. A. Akimov, Soft Matter , 2357 (2016). [47] M. M. Terzi and M. Deserno, J. Chem. Phys. , 084702 (2017). 6 [48] M. M. Terzi, M. F. Ergüder, and M. Deserno, J. Chem. Phys. , 164108 (2019). [49] L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics (Pergamon, Oxford, 1975). [50] Y. C. Fung,
A First Course in Continuum Mechanics (Englewood Cliffs, NJ: Prentice-hall, 1977). [51] M. Itskov and N. Aksel, Int. J. Solids Struct. , 3833 (2004). [52] M. Deserno, Chem. Phys. Lipids , 11 (2015). [53] L. F. Braganza and D. L. Worcester, Biochemistry , 7484 (1986). [54] S. F. Scarlata, Biophys. J. , 334 (1991). [55] M. M. Terzi, M. Deserno, and J. F. Nagle, Soft Matter , 9085 (2019). [56] M. M. Terzi and M. Deserno, in Role Mech. Study Lipid Bilayers , edited by D. J. Steigmann (Springer, Cham, 2018). [57] J.N. Reddy,
Theory and Analysis of Elastic Plates and Shells, Second Edition (CRC, Boca Raton, FL, 2006). [58] J. N. Reddy,