A consistent quadratic curvature-tilt theory for fluid lipid membranes

Lipid membranes separate living cells from their environment and, in the case of eukaryotes, subdivide them further into compartments with distinct biochemical properties. These membrane-based structures usually exhibit highly intricate morphologies, and it was indeed their distinct appearance in early light- and electron microscopy studies that led to the identification of organelles such as the endoplasmic reticulum, the Golgi apparatus, mitochondria and chloroplasts, and many types of transport-related endosomal vesicles.^1–3 Today, a prime challenge of cellular biophysics is to understand how the interplay between membrane elasticity and a host of membrane remodeling proteins gives rise to this vast richness in shape and form.⁴ 

In this article, we will be concerned with the elasticity aspect of membrane biophysics, on scales ranging from the macroscopic all the way down to the nanometer level. Specifically, we will systematically construct a theory that accounts for both geometric curvature and lipid tilt and is complete at the quadratic level. We build on the seminal work of Hamm and Kozlov (HK),⁵ as well as a recent extension thereof,⁶ but we will find that a consistent quadratic theory still requires additional terms, which either had been discovered in the past but gone AWOL or had not yet been noticed. These new terms owe their existence to a more subtle treatment of both lateral and transverse shear. We will also show that the bottom-up construction of such a theory presents us with strong constraints on its elastic parameters at the emergent geometric level, which not only helps to determine their values but also affords opportunities to test the theory.

Due to the quasi-two-dimensional structure of lipid membranes, their energy can be expressed purely in terms of their surface geometry—at least on length scales about an order of magnitude larger than membrane thickness.^7–9 The resulting curvature elastic theory, widely referred to as the “Helfrich Hamiltonian,” has been a cornerstone of lipid membrane science for more than 40 years.^10–12 It expresses the surface energy density of a membrane in terms of two geometric scalars constructed from the two local principal curvatures c₁ and c₂, namely, the total extrinsic curvature K = c₁ + c₂ and the Gaussian curvature K_G = c₁c₂,

where κ and $κ¯$ are called bending modulus and Gaussian curvature modulus, respectively, and K₀ is a membrane’s spontaneous curvature (which in most artificial systems vanishes due to up-down symmetry).

Owing to a fascinating mathematical quirk of two dimensions, the integrated Gaussian curvature is a topological invariant (up to a boundary term),^13,14 and so a closed membrane’s Gaussian curvature energy is shape independent and only changes during topological transitions, specifically by $4πκ¯$ during fission and $−4πκ¯$ during fusion. However, the magnitude of these discrete shifts is hard to pinpoint since the value of $κ¯$ is very difficult to measure—precisely because it largely hides in topology. Most studies place it at $κ¯≈−κ$⁠,^15–23 in the middle of the range [−2κ, 0] required for the quadratic expression in Eq. (1) to be positive definite in {c₁, c₂}.²⁴ 

Instead of directly measuring $κ¯$⁠, one might try to utilize a microscopic underpinning of Helfrich theory to estimate $κ¯$ in terms of some more fundamental elastic description. A well-known prediction of such a bottom-up approach is that the Gaussian curvature modulus can be obtained as the second moment of a bilayer’s lateral stress profile σ₀(z),^5,25,26

Even though the lateral stress profile cannot be measured in experiments, it is relatively straightforward to calculate in computer simulations, apart from the fact that a local definition of a stress tensor remains conceptually subtle,^27–29 and injudicious choices affect the calculation of a membrane’s lateral stress profile.^30–32 However, $κ¯$ values that are measured from the second moment generally do not match the results obtained from a patch-closure procuedure,^22,23 and they frequently even violate the stability condition $κ¯<0$⁠.^23,33–36 These discrepancies motivate us to revisit the type of bottom-up elastic theories on which Eq. (2) rests.

Among the theoretical frameworks proposed to explain the Helfrich Hamiltonian in more fundamental terms, we will focus on the influential theory by Hamm and Kozlov (HK).⁵ Its underlying starting point is itself an elastic model that permits a highly predictive analytical treatment while making very few specific assumptions: a lipid monolayer (i.e., a single bilayer leaflet) is described by the most general three-dimensional elastic energy functional that is quadratic in the strain tensor, satisfies all extant symmetries, but also couples the strain linearly to a pre-existing lateral stress profile σ₀(z). Out of this arises not only classical curvature elasticity but also a novel microscopic degree of freedom, whose coupling to the geometry is fully specified: lipid tilt—the average deviation of a lipid’s orientation away from the local membrane normal.

In a recent paper, Terzi and Deserno (TD)⁶ revisited HK theory and argued that it misses a coupling term between curvature and (the divergence of) tilt. Remarkably, the prefactor of this term is the second moment of the lateral stress profile, but taken only across a single leaflet and centered around that leaflet’s pivotal plane. Within both HK theory and TD theory, this is nothing but the Gaussian curvature modulus of a monolayer, $κ¯m$⁠. On the one hand, this is very exciting because it liberates this modulus from its purely topological provenance and suggests ways for accessing it through far more conventional means, such as the power spectrum of membrane undulations.⁶ On the other hand, it offers no answer to the aforementioned conundrum that the values of Gaussian moduli derived from second moment conditions appear dubious.

In this paper, we extend and quadratically complete HK theory and TD theory by two additional terms, both of which arise from subtleties associated with the two different ways in which shear arises in HK theory. First, lateral shear is inextricably tied to the notion of fluidity, but this notion affords a local and a global interpretation, with the latter giving rise to a new degree of freedom: lipid twist. This possibility had actually been fully embraced by Hamm and Kozlov,⁵ but it was omitted in subsequent treatments, likely because no experimental or numerical evidence existed that called for this more involved interpretation. We will see, though, that this term gets us a long way toward solving the second moment mystery of the Gaussian curvature modulus. Second, transverse shear is the physical origin of lipid tilt. Its magnitude has been treated as constant along the (very thin!) transverse direction, but we will show in the present paper that this is nevertheless not fully consistent. Permitting its value to change across the leaflet—an immediate consequence of the fact that the local normal vectors of translated lipid reference surfaces do not remain collinear either—creates yet another curvature-tilt coupling term. Fortunately, its form is identical to the one recently discovered by Terzi and Deserno,⁶ and so the amended theory has the same structure compared to TD. However, the meaning of the prefactors changes, and hence the connection between the form of a bilayer’s power spectrum and the monolayer’s Gaussian curvature modulus is more complicated than initially claimed in Ref. 6.

With the addition of these two terms, we arrive at a quadratically consistent curvature-tilt elastic theory for lipid membranes. Since it is not just postulated in the spirit of an effective field theory, whose phenomenological prefactors are not constrained beyond the exploited symmetries, we get strong constraints on its coupling constants that arise from the underlying postulated ontology. We will show that these not only provide means for estimating the values of these coupling constants but also afford opportunities to test the validity of our underlying microscopic assumptions. To illustrate both, we will reanalyze a large set of 12 different lipid membrane systems, which were extensively studied in a recent paper by Venable, Brown, and Pastor (VBP),³⁶ and whose atomistic simulation trajectories were kindly provided to us by the authors.

The large-scale bending energy of membranes, Eq. (1), features three elastic parameters: bending modulus κ, Gaussian curvature modulus $κ¯$⁠, and spontaneous curvature K₀. There are two conceptually different ways to determine these: (i) directly measure them in experiments or simulations, and (ii) start with a more microscopic theory and unearth the values of ${κ,κ¯,K0}$ in terms of the parameters of this underlying theory. These might not be known either (there lurks indeed an infinite regress), and they might not even be smaller in number; but such a procedure conceivably unveils connections between parameters at the larger scale that might not be evident at that scale per se. This is especially valuable once we refine the large-scale theory and include finer observables—such as, for instance, lipid tilt—that permit multiple novel couplings to the existing degrees of freedom. All of these require their own coupling constants, but the constraints enforced by the bottom-up approach forestall an unchecked proliferation of empirical parameters in the large-scale theory.

In their seminal work, Hamm and Kozlov⁵ proposed a bottom-up approach for the curvature-tilt elastic theory of lipid monolayers. Their derivation starts from the most general three-dimensional quadratic (in terms of, say, Green–St-Venant strains) energy density and determines the two-dimensional energy density using three main assumptions. The first one is in-plane rotational symmetry. The second one is a thin-plate assumption, which takes the transverse stress to vanish throughout the membrane and with the additional (but excellent) assumption of incompressibility helps to re-express the normal strain in terms of the lateral one. The third assumption is lateral fluidity, which sets the lateral shear modulus to zero. We will see that, as foreshadowed in HK, a subtle difference arises depending on whether this assumption is made locally or globally.

After exploiting these symmetries, the three-dimensional energy density (without the modifications to the lateral and the transverse shear we will soon discuss) is given by⁵ 

where ϵ(z) is the (depth-dependent) area strain, z is the transverse coordinate in the flat monolayer, and T_p is the tilt field on the pivotal plane—the reference surface on which the area element is independent of curvature and tilt deformations. To emphasize that the value of a position-dependent field variable is taken at that location, we will use the subscript “p.” Moreover, the superscript “0” on $e3d0$ serves as a reminder that this is the energy density used in previous studies (HK as well as TD) to which we will add extra terms in the present paper. Finally, the elastic moduli $Ẽ$(z), λ_t(z), and the lateral stress profile σ₀(z) are depth-dependent material parameters and assumed to be given. This three-dimensional energy density can be turned into a two-dimensional one by integrating over the transverse direction.

The area strain ϵ(z) arises from the deformation of the membrane and its constituent lipids. It should thus be expressible through a leaflet’s surface geometry (captured by the metric g_αβ and the curvature tensor K_αβ) and the lipid’s tangential tilt field (given by T = T^αe_α, where e_α is the coordinate basis and T^α are the tilt field’s contravarient components). In a recent publication, TD⁶ amended HK’s derivation of this connection by a term that involves curvature gradients,⁶ 

Here, $K̃$ and $K̃G$ are the effective total and Gaussian curvature, respectively, defined as the trace and determinant of the effective curvature tensor,

which leads to

where ∇_α is the metric-compatible covariant derivative and ε^αβ is the contravarient Levi-Civita tensor.^37,75 We note that another common name for $K̃$⁠, preferred by HK, is “splay.”

With the addition of a lateral shear term and a z-dependent tilt field, the energy density [Eq. (3)] generalizes to

where λ_s(z) is the depth-dependent lateral shear modulus and u_ij are the components of the three-dimensional strain tensor. Observe that the particular combination of its elements in the last term of Eq. (8) is the quadratic lateral shear strain, which is also equal to the contraction of the deviatory part of the in-plane strain tensor with itself.

There are two key differences between Eqs. (3) and (8): First, we have no longer assumed that the tilt field is constant in the transverse direction, T_p → T(z), and second, we added a lateral shear strain term [the second line of Eq. (8)]. Notice that this contribution vanishes identically if one simply assumes λ_s(z) ≡ 0, but this is not the most general way in which to enforce lateral fluidity in an effective 2d theory.

This is the gist of the differences between HK/TD and our present paper. To support this with a few more technical details, let us in Sec. II A compare the two fluidity assumptions and explore their implication for lipid elasticity. In Sec. III B, we then examine the consequences of the z-dependence of the tilt field.

Fluid materials offer no resistance to a constant shear stress, and this is incorporated into elastic theories by setting the associated shear modulus to zero.⁷⁵ In our case, membranes are modeled as 2d fluid surfaces, where the lipids are free to move (“flow”) in-plane but restricted along the out-of plane direction. However, there is a subtlety in the present situation related to the procedure by which the final effective 2d theory is constructed. Its energy density is determined after dimensional reduction, i.e., after integrating the 3d energy density (8) over the transverse z-coordinate. Applying the fluidity assumption before or after this process leads to different results because there is more than one way to make the z-integral vanish. In fact, HK suggest two ways to apply the in-plane fluidity assumption, which can be summarized as follows:

In the first case, the lateral shear modulus vanishes locally (i.e., at every z-value) in the membrane, and hence the lateral shear term drops out in its entirety from the energy functional. This is more restrictive than the second way of applying in-plane fluidity, in which only the zeroth moment of λ_s(z) vanishes (which couples to the lateral shear in the effective 2d theory). The difference is that in this second case, higher moments of the lateral shear modulus, which arise in the subsequent algebra, survive and leave their trace in the effective 2d theory. To see how this happens, let us rewrite the shear term using surface geometry and tilt field,⁵ 

Since the expressions in Eq. (10) are multiplied by λ_s(z) in Eq. (8), dimensional reduction will hence produce three additional terms in the effective two-dimensional energy functional, all proportional to the second moment of the lateral shear modulus profile. The functional form of the first two (involving $K̃2$ and $K̃G$⁠) already appears in the existing theory, and so this addition simply renormalizes the value of their respective coupling constants. But the third term is new. Since the curl of the lipid tilt field is usually called lipid twist, we see that the third term quadratically penalizes lipid twist, and this reveals the second moment of the shear modulus profile to be equal to the lipid twist modulus.

Let us provide an illustration for how a different take on membrane fluidity lies at the heart of twist rigidity. Consider the small patch of a single lipid leaflet sketched in Fig. 1. The lipids in the front row tilt to the right, while those in the back row tilt to the left, and this tilt changes continuously from front to back. If the entire lipids were to move as their heads do, this would be a horizontal shear deformation, and it would be energetically irrelevant due to fluidity. But the tails move oppositely to the heads, thereby “undoing” the shear of the heads. Nevertheless, this net-shear-free deformation amounts to a twist deformation as is obvious by looking at the colored strip of lipids in the middle of the picture. The lipid twist term penalizes such a deformation, but it can only do so if the shear modulus does not identically vanish.

Recent molecular dynamics simulations have shown that the power spectra of lipid orientation fluctuations at a large wave vector require for their quantitative description a lipid twist penalty in the energy functional.^36,38,39 This suggests that the second moment of the lateral shear profile does indeed not vanish, and the more realistic second fluidity assumption is closer to the truth—a piece of empirical information not available to Hamm and Kozlov when they developed their theory. What makes this so interesting is that this lateral shear term does not merely elevate twist to a relevant player; it also changes the definitions of both bending modulus and Gaussian curvature modulus in terms of the microscopic parameters.⁵ We will come back to this point and discuss the implications of this modification to simulation results in Sec. III.

The quadratic tilt term in Eq. (8) originates from the transverse shear component of the strain tensor. In shear deformation plate theories, this mode describes the following phenomenon: consider the underformed plate and “paint” a line through it—along its local geometric normal, but anchored to the material; as the plate is bent, these two directions will progressively deviate from one another.⁴⁰ For our case of a lipid leaflet, this finds a natural correspondence to the tilting of the lipids’ orientation away from the local membrane normal. If we write n for the normalized lipid director field and N_p for the geometric surface normal on the pivotal plane, we can define a tilt vector T_p as the difference,

where we divide the first term by n · N_p to enforce T_p · N_p = 0, so tilt is purely tangential. If θ is the angle associated with this tilt, we may easily check that tan θ = |T_p|.

Curvature-induced tilt is generally small: even at very strong bending, θ is at most a few degrees.⁴¹ Hence, the presumably small variations of an already small observable have so far been ignored.^5,6 To be sure, the reason why tilt would even change as a function of position within the membrane is very subtle. The lipid director, n, is constant along the lipid, unless we permit the lipid itself to bend—a complication we will indeed not pursue further (but see our comments below). However, the same is not true for the geometric normal N(z). The crucial point is that the material reference surfaces that foliate the membrane into layers are not merely translations of the pivotal plane along its local normal (i.e., they are not merely parallel surfaces); if they were, their normals would indeed coincide with those of their parent surface.¹⁴ Instead, our material surfaces arise from translations along the local lipid orientation n (Ref. 6 hence calls them “lipid-shifted surfaces”). Unfortunately, this renders their associated geometric normals position dependent, and so tilt should technically be defined as

where each value of z represents a different lipid shifted surface. In the  Appendix, we show that this creates a change in the tilt field away from the pivotal plane that is quadratic in the distance ζ measured along the lipid in the curved configuration and proportional to the gradient of the local curvature Eq. (A6). This, in turn, changes the square of the tilt according to

Maybe somewhat surprisingly, the correction to the bare $Tp2$ is not negligibly small: it is proportional to the curvature and the tilt field and in fact of exactly the same form as the correction term in the area strain of Eq. (4). Fortunately, this correction to the tilt field spawns no terms beyond this instance because everywhere else in the Hamiltonian tilt already refers to the pivotal plane.

Having determined the z-dependence of the tilt field, we can write it (and all expressions constructed from it) in terms of quantities defined on the pivotal plane. Hence, these are the only ones left after the dimensional reduction (i.e., after integrating over z to get the surface energy density), and so we can declutter our notation by dropping the subscript “p.” There is no other tilt field left with which we could confuse it.

One might wonder: why stop at lipid tilting? Could one include further deformations that allow lipids to deviate from an average straight shape—a deformation we might call “lipid flex”? If so, would such deformations leave terms at the quadratic order we aim for? The short answer is that such deformations could be included and that they would almost certainly leave quadratic terms (if nothing else, a “flex modulus” must arise), but that incorporating them is likely not a productive way forward. For reference, in continuum mechanics parlance, the kinematic constraint that material normals do not deviate from surface normals during thin plate bending is called Kirchhoff-Love theory, while the relaxed constraint that they can rotate is known as Mindlin-Reissner theory or “first-order shear deformation theory,” and it is typically only used for thick plates (say, once the ratio between thickness and planar extension exceeds 10%). Considering that even under extreme bending, lipid tilt rarely amounts to more than a few degrees, we expect the numerical corrections due to flexing to be even smaller. Moreover, since the flex modulus itself and potential other novel cross terms would arise as further moments of either λ_ijkl(z) or σ₀(z), it will be difficult to determine them ahead of time. (We will see below that even estimating κ_st causes significant trouble.) In light of this, we worry that opening a flex degree of freedom will essentially amount to the introduction of further parameters that are exceedingly hard to determine, render the systematic counting of “orders” even less obvious, and tarnish the robustness of a theory based on simple and testable assumptions. We hence exclude such deformations from our considerations.

We can now proceed with the dimensional reduction, in which the monolayer surface energy density—including the lipid twist from Eq. (10) and the correction to the tilt field from Eq. (13)—emerges from integrating the three-dimensional energy density over the transverse z-coordinate,

In the second term of the second line, we integrated by parts to shift the derivative from curvature to tilt, $−T⋅∇K̃→+K̃∇⋅T$⁠.⁷⁶ Recall that both $K̃$ and T refer to the pivotal plane, and all occurrences of z are now explicit: they lead us to calculate moments (of degree 0, 1, and 2) of various elastic moduli or of the stress profile. For later convenience, let us introduce the following notation for these moments:

where the integral spans the width of a flat monolayer, which, in practice, means the range over which the moduli or the stress profile are nonzero and—we recall—where z = 0 lies at the pivotal plane (not the bilayer midplane). With this notation, Eq. (14) can be re-expressed as

where numerous material parameters arise that can be expressed as combinations of moments,

They are conventionally named as follows: κ_m is the bending modulus, $κ¯m$ is the Gaussian curvature modulus, K_0,m is the spontaneous curvature, κ_t,m is the tilt modulus, and κ_tw,m is the twist modulus. The modulus κ_st,m is new in our theory, and since it couples splay $K̃$ to (the divergence of) lipid tilt, we propose to call it splay-tilt (coupling) modulus. Since all these parameters refer to the monolayer, we decorate them with an additional subscript “m.”

Observe that the twist modulus arises as the second moment of the shear modulus profile λ_s(z)—an object that would be identically zero, were it not for the more subtle notion of fluidity, as discussed in Sec. II A. Recall, in particular, Eq. (10), which showed that the surviving term comprises three expressions of different structures, only one of them directly related to lipid twist. Hence, a nonzero twist rigidity also changes the numerical value of two other elastic parameters: the bending modulus κ_m and the Gaussian curvature modulus $κ¯m$⁠. Notice that these microscopic underpinnings also have other interesting implications, for instance, $2κm+κ¯m$ is independent of our assumptions of lateral shear (since this combination is independent of λ_s) or $κm+κ¯m$ is independent of the stress profile.

Comparing Eq. (17) with the recent results of Terzi and Deserno,⁶ we find differences in four relations—three due to the more subtle treatment of lateral shear (i.e., fluidity) and one due to the consistent treatment of transverse shear (the corrections to T_p). Hamm and Kozlov had previously discovered the twist corrections,⁵ but subsequently discarded them, presumably for lack of experimental evidence that twist mattered. Besides the twist term [the last line in Eq. (16)], all other expressions occurred both in the HK⁵ functional and the TD⁶ functional, albeit with different microscopic underpinnings for the associated elastic parameters. In particular, for the earlier TD functional κ_st,m turned out to coincide with $κ¯m$⁠, which, in principle, would permit one to determine the Gaussian curvature modulus of the monolayer by observing tilt-related corrections to observables such as the power spectrum of membrane undulations.⁶ The emergence of two different corrections to these two moduli has made such an approach considerably more tricky.

The two individual leaflets of a bilayer can slide past one another, and so the bilayer energy density is the sum of the two individual monolayer energy densities. However, since we used the pivotal plane as the reference surface for a single leaflet, while for the bilayer, we would evidently like to use its midplane, we first have to perform a coordinate transformation that shifts the monolayer reference surface to what will become the bilayer’s midplane, i.e., essentially the “bottom” surface of a leaflet, where the lipid tails end. The key points of this calculation have been outlined in Ref. 6 to which we refer the reader for technical details. Here, we instead review in general terms what differences we should expect.

Obviously, this coordinate transformation involves a shift by the distance between pivotal plane and midplane, which we denote by z₀ in the flat bilayer. Any correction term involving it must therefore be at least of order z₀ and hence—for dimensional reasons—contain at least one additional factor of either the curvature or the tilt gradient. The smallest correction to any quadratic term in Eq. (14) should thus be cubic in inverse length, an order which we systematically ignore.

This expectation is correct—except for the case of the squared tilt: this is the only quadratic term whose dimension is not inverse squared length but which is actually dimensionless. Luckily, the geometry underlying a shift of the reference surface is precisely the same that arose when we discussed transverse shear and a z-dependent tilt field in Sec. II B. The main result from there, Eq. (13), translates immediately, and so we get

We alert the reader that this correction term was missing in the derivations presented in TD.⁶ While not creating a novel term in the energy density, it changes the relation between the splay-tilt modulus of the bilayer, κ_st, and its monolayer partner κ_st,m.

To complete the derivation, it is once again useful to combine the two tilt-fields T_↑ and T_↓ in the two leaflets into their symmetric and antisymmetric combinations since they qualitatively differ in how they couple to the geometry,

Notice that this only makes sense after the reference point shift, for vectors living on different reference surfaces cannot be added. With this, the bilayer energy density becomes

where $K̃G−$ is the determinant of $K̃αβ−≔Kαβ+∇αTβ−$⁠. The elastic parameters of the bilayer are simple combinations of their monolayer counterparts,

The fourth expression is new; all the others agree with previous works.^25,26,38 Observe that T₋ couples directly to the curvatures (both mean and Gaussian), while T₊ does not—other than “minimally,” meaning, other than through the metric-dependence of the curvilinear divergence and curl. This is the main reason why the (anti-)symmetrization of Eq. (19) is useful.

In Sec. II, we have derived the complete quadratic curvature-tilt theory for both monolayer and bilayer membranes as it emerges from the most general quadratic elastic functional of a fluid thin film. We have, in particular, drawn attention to the subtle role played by the precise meaning of “lateral fluidity” and to the consistent treatment of the transverse shear term. Besides giving rise to lipid twist, the resulting corrections only altered the definitions of some of the membrane elastic parameters in terms of their underlying microscopic ones—see Eqs. (17) and (21).

We will now look at some of the consequences of these modifications, the most visible of which is the occurrence of new moduli, κ_st,m and κ_st. In TD theory, these coincided with the Gaussian moduli and hence did not require a separate name. The two additional but different terms in Eqs. (17b) and (17e) break this accidental (indeed: erroneous) symmetry.

In order to test the predictions of our new theory, we use it to analyze the power spectra of 12 different atomistic bilayer trajectories, which were kindly provided by Rick Venable and Richard Pastor and which had been previously analyzed by Venable, Brown, and Pastor (henceforth VBP) in Ref. 36. Before delving into a discussion of the results, let us first, for completeness, revisit the expressions following from our theory with which these power spectra were fitted. Their derivation follows with only minimal modifications from Ref. 6 and will not be repeated here.

Among the elastic parameters for lipid membranes, the bending modulus κ is arguably the most important one and probably the best studied. It is typically obtained—both in simulations^19,42–57 and in experiments^58–64—by analyzing the power spectrum of thermal shape undulations, i.e., the mean squared amplitude ⟨|h_q|²⟩ of a membrane’s Fourier modes h_q. However, Watson et al.⁶⁵ showed that one can also obtain κ by analyzing the longitudinal component of the director fluctuation spectrum, $⟨|nq∥|2⟩$⁠, and that this analysis works well for membrane sizes smaller than the ones typically required for the height fluctuations. Adapting to the notation used in the present paper, the splay-tilt coupling term introduces a divergence at a critical wave vector q_c to the longitudinal component of director fluctuations,⁶ 

where

and

On the other hand, the tilt and the twist moduli can be obtained from the transverse component of the director fluctuation spectrum,⁶⁵ 

Notice, in particular, that the transverse director fluctuation spectrum becomes flat (i.e., q-independent) if the twist modulus vanishes. This clearly contradicts the measured spectra (see the supplementary material) and hence supports the more subtle global fluidity assumption (9b) on which the existence of a nonzero κ_tw relies.

The power spectra from the director fluctuations are hence informative about four elastic parameters: {κ, κ_st, κ_t, κ_tw}. Observe that due to the additional corrections discussed in the present paper, this set contains κ_st and not $κ¯m$⁠, as in TD.

VBP analyze atomistic simulations (based on the CHARMM36 force field) of 12 different one-component lipid membranes.³⁶ Table II lists a few of the observables that were extracted from measurements unrelated to fluctuation spectra; we will use them in our discussion below. Here, we reanalyze the simulation trajectories of VBP and fit them to the spectra in Eqs. (22) and (23). More specifically, we simultaneously fit the longitudinal and transverse components of director fluctuations to Eqs. (22a) and (23) in order to obtain a set of {κ, κ_st, κ_t, κ_tw} values for each of the 12 different lipid types. Details of this procedure, as well as graphs of all power spectra, and their associated fits can be found in the supplementary material. Due to the way in which κ_st enters the fluctuation spectrum (22a), both signs of κ_st result in exactly the same spectrum and hence cannot be distinguished. Since in Sec. III B 2 we will derive estimates for κ_st using the insight that we have from bottom-up theory, which for all but one lipid type (namely: palmitoyl sphingomyelin) we find to be positive, we choose the positive sign for all 12 κ_st values. Results of our analysis are given in Table I; errors are determined by parametrically bootstrapping the spectra.^66,67

We have previously shown that the undulation spectra ⟨|h_q|²⟩ are better captured by our theory, which includes the high-q divergence, as compared to HK theory, which lacks the 1/Δ_q prefactor (even if we arrange for the same number of fitting parameters).⁶ This comparison relied on very accurate CG data, for which also the low-q part of the spectrum was present and well equilibrated. In contrast, the VBP data we use in the present paper were obtained for very small systems (L ≈ 14.4 ± 0.5 nm), which were never meant to be analyzed for their undulation spectra. We hence cannot reasonably use them to assess the relative merits of our theory with regard to its ability to describe ⟨|h_q|²⟩. As far as the lipid orientation spectra are concerned, we remind the reader that the transverse orientation fluctuations, $⟨|nq⊥|2⟩$⁠, are described by the same formula in our present theory and in HK theory [assuming HK is also amended by adopting the more subtle fluidity assumption from Eq. (9b)]. The longitudinal orientation fluctuations, on the other hand, differ: our theory again includes the divergence 1/Δ_q, which is missing in HK theory. This implies that $q2⟨|nq∥|2⟩$ is a constant within HK theory, while the actual data show it to increase for larger q—a feature captured by our theory, which is indeed informative about q_c, and hence the combination $κst2/κκt$⁠. An even more detailed comparison would require us to also look into systematics (e.g., questions how the results depend on the precise definitions of lipid directors or the number of q-vectors we fit to), which we defer to a future publication.

For completeness, let us conclude this section with a brief comment on the undulation spectra. We have seen that fitting the director fluctuations results in two values for κ_st (and hence r_st), which differ in their sign, so that only their modulus can be determined. A similar degeneracy arises in the height fluctuation spectrum,⁶ but with an interesting difference. Recall that

where

Again, two sets of {κ_t, κ_st} values exist that result in exactly the same spectrum, so there is no way to differentiate between them on the basis of ⟨|h_q|²⟩ alone. The reason is that the transformation from the immediate fitting parameters {κ, κ_t,eff, q_c} to the derived set {κ, κ_t, r_st} involves solving a quadratic equation for r_st,

with

Obviously, r_st+ > 0 and r_st− < 0. Moreover, r_st− ≥ −1, while r_st+ ≤ 1 only when Q² ≥ 8. Notice also that the two r_st values are not just negatives of each other, for they genuinely differ in their magnitude: r_st+ > −r_st−. Comparing this with the degenerate r_st values obtained from the analysis of transverse orientation fluctuations implies that one of the two overall signs must agree better and might hence help to break the degeneracy. Unfortunately, given that our present data do not allow a reliable undulation analysis, this option is presently not available to us.

As a final comment, Wang and Deserno have previously shown that the tilt modulus can also be extracted from simulations of buckled membranes.⁴¹ Since their analysis relied on the original HK Hamiltonian, TD subsequently pointed out that the tilt modulus $κtWD$ extracted following the Wang-Deserno prescription must be multiplied by the correction factor $(1+12rm)$⁠, where $rm=κ¯m/κm$ is the ratio between a monolayer’s Gaussian and mean bending modulus.⁶ This factor arises as a consequence of the new coupling term discovered in TD, but in the present paper, we argue that this coupling term is not multiplied by the prefactor $κ¯m$ but by the new modulus κ_st. This, however, merely changes the correction factor to $(1+12rst)$⁠. Observe that the occurrence of lipid twist—ignored in both Refs. 41 and 6—causes no further modifications because the essentially one-dimensional deformation of a buckle does not excite lipid twist, and hence twist drops out of the Euler-Lagrange equation of a buckled membrane.

In Sec. II, we expressed the elastic constants of the curvature-tilt theory in terms of the underlying material properties of a prestressed thin fluid elastic film. The monolayer elastic constants are given in Eq. (17), and the bilayer counterparts are given in Eq. (21). It is now time to look more closely at these expressions and explore several relationships they imply between the elastic constants.

Let us begin with the bending modulus κ, which in terms of the underlying mechanical parameters is given by

a relation that follows by combining Eqs. (17a) and (21a). Recall that the three terms on the right are the second moments of the z-dependent Young modulus, lateral prestress, and lateral shear modulus, respectively, taken over a single leaflet and centered at the pivotal plane. Both the stress- and the shear moment can be extracted from simulations because the stress profile σ₀(z) is measurable,^30,68,69 and the shear moment is proportional to the twist modulus, which can be extracted from lipid orientation fluctuations [see also Eq. (23)].^5,38 However, there is no direct way to extract $[[Ẽ]]2$⁠. To make headway, we will assume that the Young modulus profile is actually a constant, $Ẽ$(z) = $Ẽ$⁠, which implies

where d_m is the monolayer thickness.

From an atomistic perspective, the boundaries of the integral in Eq. (27a) are not particularly well defined.⁷⁰ Here, as well as in all other moment integrals, they account for the fact that the observables we integrate over drop to zero outside a leaflet, in a way that [except for σ₀(z)] we generally do not know. Assuming that, say, the Young modulus is not z-dependent, while also imposing hard boundaries, amounts to approximating $Ẽ$(z) by a box function, which then comes down to a particular choice for a Gibbs dividing surface between bilayer and water. While the lower boundary obviously coincides with the bilayer midplane, which can be defined without much trouble, the upper boundary might, for instance, be chosen as the Luzzati plane,⁷¹ which replaces the smooth lipid/water partitioning by one sharp volume-preserving divide. Still, the choice is not unique, and so it is appropriate to treat d_m as an effective thickness—essentially, a model parameter.

Assuming volume incompressibility, a membrane’s area expansion modulus is the zeroth moment of its Young modulus,^24,72 and so the assumption of a constant $Ẽ$ can be further used to calculate

Using this, we can eliminate $Ẽ$ from Eq. (27b) and, together with Eq. (26), arrive at an estimate for κ in which each term is, in principle, accessible,

To check how well this estimate fares, we calculate bending rigidities for 12 different lipid membrane systems, using the parameters calculated by VBP (see Table II),³⁶ and compare with bending rigidities extracted from the same simulation trajectories (kindly provided by VBP), using our slightly amended mode spectrum theory from Sec. III A. The power spectra and the fits can be found in the supplementary material.

VBP provide two different monolayer thicknesses: half the trans-bilayer separation between phosphate groups, d_P:P, a proxy of the Luzzati width, and half the hydrophobic thickness, taken to be the distance d_C2:C2 between the C₂ carbons of lipid’s chains. As we discussed earlier, the intrinsic ambiguity of the notion of “membrane thickness” leads us to define an effective thickness,

We require the scale factor f^est to be the same for all 12 different lipids.

Table III lists κ values resulting from reanalyzing the VBP trajectories, together with predicted rigidities κ^est calculated via Eq. (29) and their relative difference. We determined the effective thickness coefficient by minimizing the mean square difference between measured κ and estimated κ^est values, finding f^est ≈ 1.13. Figure 2 illustrates the agreement between estimated values and measured values. The correlation coefficient between κ^est and κ is found to be r = 83%.

We will now compare the estimate of κ in Eq. (29) with two other predictions, namely, the polymer brush model⁷³ and the uncoupled two-dimensional fluid film model.⁴² The polymer brush model predicts a simple relationship between bending modulus and area expansion modulus,⁷³ 

where d_h is the thickness of the hydrophobic region. For the latter, we could use d_C2:C2, but recalling the slight arbitrariness of thickness definitions and to also make a comparison more fair, we will again define an effective distance $dheff=fPBMdC2:C2$⁠, where the scale factor f^PBM is again determined by minimizing the squared deviation between κ^PBM and κ over the 12 different lipids. This leads to f^PBM ≈ 1.26, showing that the “bare” hydrophobic distance slightly underestimates the best choice of membrane thickness. The associated rigidities κ^PBM and their relative difference from κ are also listed in Table III, and a correlation plot is shown in Fig. 2; the correlation coefficient, r = 63%, is significantly smaller than for κ^est.

A third relation between κ and K_A, derived by Goetz, Gompper, and Lipowsky,⁴² results if one views bilayers as two uncoupled two-dimensional fluid sheets, leading to

where again d is some version of the bilayer thickness. If we again introduce a scale factor f^GGL and use d_C2:C2 as the base reference, we immediately see that this is equivalent to the polymer brush model since the factor-of-2 difference in the denominator can be absorbed into f^GGL. The conceptual details of this model make it more plausible, though, that d should be identified with the Luzzati distance, and so we will use d_P:P as the base and define d = f^GGLd_P:P. Minimizing the squared deviation between κ^GGL and κ across all lipids then leads to f^GGL = 1.32. The results together with the relative differences with respect to κ are again shown in Table III and Fig. 2. The correlation coefficient is r = 58%. This differs from the polymer brush correlation because the ratio d_P:P/d_C2:C2 is lipid-dependent.

Finally, the value of f^GGL shows that its effective d is even bigger than the Luzzati width of a bilayer or, alternatively, that the denominator 48 is too large. A better value would be $48/(fGGL)2≈40$⁠. Interestingly, a prediction for κ resting on the sliding of two thin elastic sheets would give²⁴ 

where ν is the Poisson ratio of the material, which has the value $ν=12$ in the incompressible limit. Hence, this formula requires d to be about 5% smaller than the Luzzati width, while the GGL formula requires d to be about 10% bigger. The correlation coefficient is of course the same for both cases, 58%, since both expressions refer to the same base width, d_C2:C2.

Overall, we find that the estimates based on Eq. (29) correlate significantly better with the measurements of κ than those given by the polymer brush theory, the two-dimensional fluid sheet model, and the two sliding elastic thin sheets model, even if in all cases we permit a scale factor to optimize the notion of membrane thickness, which should be expected to be different in all models. The better agreement comes at the prize of a significantly more complex formula, requiring input beyond basic elastic and structural properties. This, of course, is also a lesson in the underlying elastic origins of lipid bilayer bending; Eqs. (26) and (29) show that the bending modulus reflects three distinct pieces of physics: membrane stretching, differential work against prestress, and lipid twist. Since for different lipids these three terms are not necessarily all just proportional to one another, any prediction that merely relies on one of the three (here: K_A) is likely to miss subtle variations between lipids that differ in their architecture.

Even though the predictions of Eq. (29) are overall fairly good and noticeably better than those of the (admittedly simpler) predictions from Eqs. (31) and (32), a few strong outliers remain. For instance, the bending modulus of POPG is severely underpredicted, while that of SDPE is strongly overpredicted. Is it possible to attribute these deviations to specific complications with a given lipid? As it turns out, this is very difficult because the lipids in our set were carefully chosen by VBP to cover a broad set of conditions, and so almost every lipid has a special property that could lie at the heart of some specific problem. For instance, POPG is the only charged lipid. In other cases, seemingly minor changes appear to have a big effect; for instance, SDPE and PDPE differ only by two carbons in their saturated chain, and yet the rigidity of SDPE is significantly overpredicted, while that of PDPE is somewhat underpredicted. At this point, no obvious pattern strikes us; and at any rate, it would then have to be independently verified by additional simulations of lipids not belonging to the test set of lipids from which any such pattern was extracted.

The coupling modulus κ_st enters the longitudinal director fluctuation spectra, $⟨|nq∥|2⟩$⁠, and we have used this fact to obtain κ_st by simultaneously fitting them to the analytical expressions in Eqs. (22a) and (23). Also, Eqs. (17e) and (21d) give a definition of κ_st in terms of the two-dimensional elastic parameters,

where z₀ is the distance between the midplane and the pivotal plane. In order to estimate κ_st, we reprise our idea from Sec. III B 1 and assume a constant modulus profile—this time for the tilt modulus λ_t(z) = λ_t; this yields

where d_m is the monolayer thickness. In this case, the tilt modulus becomes

An estimate for the coupling modulus can now be written in terms of the second moment of the stress profile, the tilt modulus, the monolayer thickness, and the position of the pivotal plane in the monolayer,

The values of κ_st obtained from fitting the fluctuation spectra and those from the estimate (37) are listed in Table IV. Unfortunately, the agreement is not as good as in the case for the bending modulus κ. The correlation coefficient between κ_st and $κstest$ is found to be −0.11, which means the prediction is unsuccessful. To be sure, it gets the range of magnitudes approximately right, but this is a low bar for comparison.

It is not immediately clear why the prediction for κ_st fails, while the one for κ works very well. Of note, κ refers to a low-q observable and can be determined from the low-q limit of the tilt spectra, in which microscopic confounding factors must vanish. The same is not true for κ_st, which is a high-q observable and is obtained from the high-q region of the fluctuation spectra, where it is much harder to disentangle from other complications of microphysics. (If we look at the length scale corresponding to the spectral divergence, 2π/q_c ∼ 2.6 nm, we realize this is about 3 lipid distances or one bilayer thickness.) Could it be that the assumption that λ_t(z) does not depend on z fails? Yes, but then why does the assumption that $Ẽ$(z) does not depend on z work? To answer these questions would require us to understand the z-dependence of various components of λ_ijkl, which in the absence of a microscopic theory of lipid bilayers is not an option.

The Gauss-Bonnet theorem states that the surface integral over the Gaussian curvature K_G can be written as a topological constant plus a boundary term.^13,14 One may hence expect that $κ¯m$ only affects the elastic energy of a membrane if it changes either its topology or its boundary (assuming it even has one). This renders measurements of the Gaussian curvature modulus very tricky.

For a lamellar phase to be stable, the bilayer elastic ratio must be in the range $κ¯/κ∈[−2,0]$⁠, for otherwise $12κK2+κ¯KG$ is not bounded below.²⁴ Even though several studies give reasonable $κ¯/κ$ values, surprisingly many simulation results are positive (and hence outside this permissible range) or are at least unexpectedly close to zero.^{22,23,33–36} Moreover, Hu et al.^22,23 proposed the so-called patch-closure method as a way to measure $κ¯$ in simulations. It exploits the fact that the closure probability of open patches involves predictable changes in a membrane’s boundary, and hence it reasons about the Gaussian curvature modulus at the level of the Helfrich Hamiltonian, i.e., without making microscopic assumptions about its elastic parameters. These authors then compared their findings with the stress profile formula (2), the “classical” microscopic underpinning and observed that the results of these two methods do not agree.

In Sec. II, we showed that the choice of the fluidity assumption changes the definition of the monolayer Gaussian curvature modulus in terms of microscopic elastic parameters. Assuming local fluidity, $κ¯m$ is simply the second moment of the lateral stress profile, but the more subtle global fluidity assumption amends this by the lipid twist modulus. Combining Eqs. (17b) and (17f), one finds

We stress that HK actually derived this equation, but later ignored the twist contribution.⁵ Notice that κ_tw,m > 0, and so the value of $κ¯m$ is always smaller than what would be expected based on the second moment formula alone. This helps shift the predicted values of the Gaussian curvature modulus (which tended to be too large) into the right direction.

To see how this formula fares, we will take the $[[σ0]]2$ results as listed in VBP³⁶ (which in the theoretical framework employed by that paper are identical to the Gaussian curvature modulus), and we calculate the twist modulus from the VBP trajectories, using Eq. (23). Table V lists the results, which are highly encouraging: for all lipids, the values of $κ¯m/κm$ are now negative and scatter centrally within the permissible range. Of course, not having any alternative prediction to compare to, we cannot ascertain by independent means whether these numbers are correct; but not having almost half of them violate an important stability condition is definitely progress.

In this paper, we have revisited a microscopic approach for justifying elastic lipid membrane theories that include both surface curvature and lipid tilt. We started in the footsteps of Hamm and Kozlov,⁵ with a general quadratic theory of a laterally fluid and prestressed three-dimensional (but thin) elastic sheet, and derived a complete quadratic curvature-tilt theory. Underscoring the importance of enforcing lateral fluidity only at the level of the resulting 2d theory, not locally (i.e., at all depths) for the 3d thin sheet, we get three additions to the curvature energy density: two which renormalize the bending and Gaussian curvature modulus and one which generates a new term that penalizes lipid twist. This, in particular, implies that the microscopic definition of the (monolayer) Gaussian curvature modulus is not merely the second moment of the (monolayer) stress profile; one also needs to subtract twice the monolayer’s twist modulus. All this was well known to Hamm and Kozlov,⁵ but in the early aughts, no evidence for lipid twist existed, and so the more involved fluidity assumption appeared to have been fallen from grace. Yet, it offers a natural solution to the (by now) well-known puzzle that Gaussian curvature moduli calculated exclusively via the second moment of the stress profile tend to be too large—often so much so that they violate elementary stability conditions. The twist modulus correction, which always subtracts a positive value, goes a long way toward fixing this problem.

We also discover a second correction, not present in earlier curvature-tilt theories:^5,6 a quadratic term proportional to both tilt and curvature gradient (or, equivalently, curvature and tilt divergence), which arises from a more careful analysis of the transverse shear strain. This term originates from the depth-dependence of the transverse shear and has the same functional form as the curvature-tilt coupling found recently by Terzi and Deserno.⁶ As such, it leaves the form of their macroscopic theory intact, but it does change the “microscopic makeup” of the coefficient multiplying the novel $K̃∇⋅T$ term. Regrettably, this prefactor used to be the monolayer Gaussian curvature modulus, but in the new version, it gets amended by the second moment of the tilt modulus profile so that a direct access to the Gaussian modulus is no longer possible.

We have also determined the bilayer version of this tilt-curvature Hamiltonian. This required shifting the reference surface from the pivotal plane to the midplane before adding the energies of two individual monolayers. For the same mathematical reasons that generated the novel transverse shear term, this reference shift also introduces a previously overlooked correction, which again changes the microscopic definition of the bilayer’s curvature tilt coupling modulus, while leaving the form of the overall theory untouched.

To test our theory against real data, we have reanalyzed the simulation trajectories studied in the paper by Venable et al.³⁶ and fit the fluctuation spectra to our predictions. A simultaneous fit of longitudinal and transverse components of director spectra reveals bending, tilt modulus and lipid twist, κ, κ_t, and κ_tw, as well as the newly introduced coupling modulus κ_st. Exploiting the fact that these moduli can be expressed in terms of microscopic parameters, we produced estimates for both κ and κ_st, which we subsequently compared to other predictions. For κ, we looked at polymer brush theory⁷³ and the two-dimensional fluid film model.⁴² Our own estimates correlate significantly better with κ-values obtained from the fluctuation spectra of 12 different lipid types than the other two estimates. However, our attempts to microscopically predict the coupling modulus κ_st were unsuccessful. This is either a sign that our theory is still missing some important physics or simply a reminder that the physics underlying κ_st happens at length scales that reach down to the molecular level, where quantitative fidelity might be too much to expect from a continuum treatment.

We provide additional details for how the director fields are defined and put on a grid, how the error in the measured power spectra is estimated, how the simultaneous fits are performed, and how their quality is assessed in our supplementary material. This file also contains graphs of all longitudinal and transverse director fluctuation power spectra, together with their fits.

We gratefully acknowledge many stimulating discussions with John Nagle, and we are thankful to Rich Pastor and Rick Venable for making the trajectories available we analyzed in this paper. M.F.E. greatly appreciates statistics advice given by Max G’Sell. M.D. also thanks the NSF for financial support via Grant No. CHE 1764257.

In this short appendix, we determine the transverse coordinate dependence of the tilt field T′ ≡ T(ζ). By definition, its components on any lipid shifted surface are given by its projection onto the tangential coordinate basis $eα′$ of that surface,

where ζ is the transverse distance measured along the lipid direction. The last term in Eq. (A1) vanishes because the normal and the tangent vectors of any surface are perpendicular. Appendix A in Ref. 6 shows that the tangent vector of lipid shifted surfaces is given by

The projection of $eα′$ onto the director field then becomes

where we used the fact that

The latter follows because n ∥ N′ + T′ and n · n = 1. The latter identity also implies n · ∇_αn = 0. Moreover, since trivially N′ · T′ = 0, multiplying Eq. (A4) by N′ leads to

This shows that the prefactor 1/(n · N′) in Eq. (A1) differs from unity only by a term that is already quadratic in the tilt field; hence, multiplying it into the numerator expression worked out in Eq. (A3) produces additional correction terms that are at least cubic in the fields and can hence be neglected. Putting everything together then gives the (quadratic order) tilt field on a lipid shifted surface in terms of the tilt field and curvature on the pivotal plane,

The quadratic tilt term in Eq. (8) is the square of this expression, T² = g′^αβT_αT_β, but this unfortunately appears to involve the inverse metric on the lipid shifted surface. Luckily, though, calculating g′^αβ turns out to be unnecessary at the order we aim for: observe that up to the quadratic order in ζ, the dyad T_αT_β is given by

This involves two distinct terms, both quadratic in the fields. Observe now that contracting Eq. (A7) with g′^αβ or with g^αβ leads to the same result (at the desired order) because the difference between g′^αβ and g^αβ is at least the first order in ζ and must hence (for dimensional reasons) involve a factor of either ∇_aT or K, i.e., it must at the lowest order be linear in the fields. Since both terms in Eq. (A7) are already quadratic in the fields, no additional terms can arise. Up to the quadratic order, we hence get

Here, we also used the relation $ζ2=z2+O(z3)$⁠. Note that the last term in Eq. (A8) has the same form as the novel term in the area strain (4). This is no coincidence: the origin of both terms is the fact that the normal vector of the lipid shifted surfaces depends on z.