Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Bilayer Elasticity at the Nanoscale: The Need for New Terms

  • Anne-Florence Bitbol,

    Current address: Lewis-Sigler Institute for Integrative Genomics, Princeton University, Princeton, New Jersey, United States of America

    Affiliation Laboratoire Matière et Systèmes Complexes (MSC), Université Paris Diderot, Paris 7, Sorbonne Paris Cité, CNRS UMR 7057, Paris, France

  • Doru Constantin,

    Affiliation Laboratoire de Physique des Solides, Université Paris-Sud, Paris 11, CNRS UMR 8502, Orsay, France

  • Jean-Baptiste Fournier

    jean-baptiste.fournier@univ-paris-diderot.fr

    Affiliation Laboratoire Matière et Systèmes Complexes (MSC), Université Paris Diderot, Paris 7, Sorbonne Paris Cité, CNRS UMR 7057, Paris, France

Abstract

Continuum elastic models that account for membrane thickness variations are especially useful in the description of nanoscale deformations due to the presence of membrane proteins with hydrophobic mismatch. We show that terms involving the gradient and the Laplacian of the area per lipid are significant and must be retained in the effective Hamiltonian of the membrane. We reanalyze recent numerical data, as well as experimental data on gramicidin channels, in light of our model. This analysis yields consistent results for the term stemming from the gradient of the area per molecule. The order of magnitude we find for the associated amplitude, namely 13–60 mN/m, is in good agreement with the 25 mN/m contribution of the interfacial tension between water and the hydrophobic part of the membrane. The presence of this term explains a systematic variation in previously published numerical data.

Introduction

As basic constituents of cell membranes, lipid bilayers [1] play an important role in biological processes, not as a passive background, but rather as a medium that responds to and influences, albeit in a subtle way, the behavior of other membrane components, such as membrane proteins [2]. The coupling between the lipid bilayer and guest molecules does not occur by the formation of chemical bonds, but rather by a deformation of the membrane in its entirety. To describe it, one must resort to concepts developed in soft matter physics for the understanding of self-assembled systems.

At length scales much larger than their thickness, the elasticity of lipid bilayers is well described by the Helfrich model [3]. However, nanometer-sized inclusions, such as membrane proteins, deform the membrane over smaller length scales. In particular, some transmembrane proteins have a hydrophobic part with a thickness slightly different from that of the hydrophobic part of the membrane. Due to this hydrophobic mismatch, the hydrophobic core of the membrane locally deforms [4][6]. As this deformation affects the thickness of the membrane, and as its characteristic amplitude and decay length are both of a few nanometers [7], it cannot be described using the Helfrich model. In fact, since the range of such deformations is of the same order as membrane thickness, one can wonder to what extent continuum elastic models in general still apply, and what level of complexity is required for an accurate description. In particular, which terms must be retained in a deformation expansion of the effective Hamiltonian?

Experimental data is available for the gramicidin channel [8], a transmembrane protein formed by two protein monomers. The channel being large enough for the passage of monovalent cations, conductivity measurements [9] can detect its formation and lifetime, which are directly influenced by membrane properties. The gramicidin channel can therefore act as a local probe for bilayer elasticity on sub-nanometer scales (see, e.g., Ref. [10]). Motivated by this opportunity, sustained theoretical investigations have been conducted in order to construct a model describing membrane thickness deformations [7], [11][13]. Recently, detailed numerical simulations have been performed, giving access both to the material constants involved in elastic models and to the membrane shape close to a mismatched protein [14][16]. This numerical data provides a good test for theoretical models.

In this article, we put forward a modification to the models describing membrane thickness deformations. We argue that contributions involving the gradient (and the Laplacian) of the area per lipid should be accounted for in the effective Hamiltonian per lipid from which the effective Hamiltonian of the bilayer is constructed, following the approach of Refs. [12], [13]. We show that these new terms cannot be neglected, as they contribute to important terms in the bilayer effective Hamiltonian. We discuss the differences between our model and the existing ones. We compare the predictions of our model with numerical data giving the profile of membrane thickness close to a mismatched protein [14][16], and with experimental data on gramicidin lifetime [17] and formation rate [18].

Results: Membrane Model

We consider a bilayer membrane constituted of two identical monolayers, labeled by and , in contact with a reservoir of lipids with chemical potential . We write the effective Hamiltonian per molecule in monolayer as(1)where is the area per lipid, while is the local mean curvature of the monolayer, and is its local Gaussian curvature (denoting by and the local principal curvatures [19] of the monolayer, we have and ). All these quantities are defined on the hydrophilic-hydrophobic interface of each monolayer. Eq. 1 corresponds to an expansion of for small deformations around the equilibrium state where the membrane is flat and where each lipid has its equilibrium area . Any constant term in the free energy per lipid is included in a redefinition of the chemical potential . From now on, we will consider small deformations of an infinite flat membrane and we will work in the Monge gauge, so and , where represents the height of the hydrophilic-hydrophobic interface of each monolayer with respect to a reference plane . The upper monolayer is labeled by and the lower one by . Many constants involved in Eq. 1 can be related to the constitutive constants of a monolayer: is the compressibility modulus of the monolayer, is its bending rigidity, is its Gaussian bending rigidity, is its spontaneous (total) curvature, and is the modification of the spontaneous (total) curvature due to area variations (see Methods, Sec. 1.1).

In the case where , Eq. 1 is equivalent to the model of Ref. [19], which is the basis of that developed in Refs. [12][15]. To our knowledge, existing membrane models including the area per lipid (or, equivalently, the two-dimensional lipid density) do not explicitly feature terms in the gradient, or Laplacian, of this variable [20]. The possibility of an independent term proportional to the squared thickness gradient was however suggested on symmetry grounds in Ref. [21], while pointing that it could arise from the specific cost of modulating the area per lipid (see note (18) in Ref. [21]). In the present work, we show that the terms in , and cannot be neglected with respect to others. We focus on the influence of , for which we provide a physical interpretation, and we will set in the body of this article in order to simplify our discussion and to avoid adding unknown parameters. However, the derivation of the membrane effective Hamiltonian is presented in Secs. 1.1–1.2 of our Methods part, in the general case where , and are all included.

The effective Hamiltonian of a bilayer membrane patch with projected area at chemical potential can be derived from Eq. 1. For this, the effective Hamiltonians per unit projected area of the two monolayers are summed, taking into account the constraint that there is no space between the two monolayers of the bilayer, and assuming that the hydrophobic chains of the lipids are incompressible. This derivation is carried out in Sec. 1.1 of our Methods part. It results in an effective Hamiltonian of the bilayer membrane that depends on three variables: the average shape of the bilayer, the sum of the excess hydrophobic thicknesses of the two monolayers, each being measured along the normal to the monolayer hydrophilic-hydrophobic interface (see Fig. 1 and Eqs. 2629), and the difference between the monolayer excess hydrophobic thicknesses. (The excess hydrophobic thickness of a monolayer is defined as the hydrophobic thicknesses of this monolayer minus its equilibrium value.)

thumbnail
Figure 1. Definitions.

A) Cut of a bilayer membrane. The solid black lines mark the boundaries of the hydrophobic part of the membrane, and the exterior, which is shaded in blue, corresponds to the hydrophilic lipid heads and the water surrounding the membrane. The hydrophobic thickness, defined along the normal to the hydrophobic-hydrophilic interface, of the upper (resp. lower) monolayer, shaded in orange (resp. yellow), is (resp. ). The height of monolayer along is denoted by . The average membrane shape, , is represented as a red dashed line. B) Cut of a bilayer membrane (with hydrophobic part shaded in yellow) containing a protein with a hydrophobic mismatch (orange square). The equilibrium hydrophobic thickness of the bilayer is , while the hydrophobic thickness of the protein is . The average shape of the membrane is flat, and the thickness deformations of the two monolayers are identical (). Hence, the average shape is constant, and confounded with the midlayer of the membrane. Although is defined along the normal to the monolayer hydrophilic-hydrophobic interface, the boundary condition at the inclusion edge, i.e., in , simply reads to first order (see main text, Section entitled “Deformation profiles close to a mismatched protein”).

https://doi.org/10.1371/journal.pone.0048306.g001

In the present work, we are not interested in the degree of freedom , which is not excited in the equilibrium shape of a membrane containing up-down symmetric mismatched proteins (see see Fig. 1B). Hence, in Sec. 1.2 of our Methods part, we integrate out, which amounts to minimizing with respect to since our theory is Gaussian. The resulting effective Hamiltonian, which involves and , is given by Eq. 32 in Sec. 1.2 of our Methods part. In this effective Hamiltonian, the variables and are decoupled, and the part depending on corresponds to the Helfrich Hamiltonian [3]. Hence, our model gives back the Helfrich Hamiltonian if the state of the membrane is described only by its average shape (see Methods, Sec. 1.3).

Here, we focus on variations of the membrane thickness, i.e., on the variable . We thus restrict ourselves to the case where the average shape of the membrane is flat (see Fig. 1B). In this case, we obtain, from Eq. 32:(2)In the case where , on which the body of this article focuses, the various constants introduced in Eq. 2 read:(3)(4)(5)(6)(7)In these equations, denotes the equilibrium hydrophobic thickness of the bilayer membrane, plays the part of an externally applied tension (see Methods, Sec. 2), is the compressibility modulus of the membrane, is its Gaussian bending rigidity, is the bending rigidity of a symmetric membrane such that , is the spontaneous (total) curvature of a monolayer, and is the modification of this spontaneous curvature due to area variations. In addition, we have introduced , which has the dimension of a surface tension, like . Note that the last three terms in Eq. 2 are boundary terms.

In Sec. 1.2 of our Methods part, the expressions of , , and are provided in the more general case where and are included.

We wish to describe a membrane with an equilibrium state that corresponds to a homogeneous thickness. A linear stability analysis (presented in Sec. 1.4 of our Methods part) shows that the flat shape is stable if , , and(8)

Discussion

Comparison with existing models

Our model Eq. 2 has a form similar to that of the models developed in Refs. [12][15]. However, it differs from these previous models on several points. First, our definition of is slightly different. Second, we have included the effect of an applied tension . Finally, the various constants in Eq. 2 have different interpretations, and thus different values, from the ones in the existing models. Let us discuss these points in more detail.

On the definition of .

In the present work, the variable , which is the relevant one to study membrane thickness deformations, is defined as the sum of the excess hydrophobic thicknesses of the two monolayers, each being measured along the normal to the monolayer hydrophilic-hydrophobic interface (see Eqs. 2629 in the Methods section). This definition of has the advantage of being independent of deformations of the average membrane shape .

The excess thickness variable used in Refs. [7], [12][15], [18], [22], [23] reads in our notations:(9)Using Eqs. 9 and 25, and working to second order, we obtain(10)which shows that there is a second-order difference between and our variable . Consequently, the difference between the definition used in the previous works and ours regards only the term linear in , i.e., the tension term, which was not included in these works. At zero applied tension, the two definitions are equivalent, i.e., it is equivalent to use or . Our definition of is the right one for rigorously taking tension into account, because it is independent of deformations of the average membrane shape : the energy stored in the variable only comes from thickness variations. (The variable of Refs. [7], [12][15], [18], [22], [23] corresponds to the difference between the bilayer hydrophobic thickness projected along and the non-projected equilibrium hydrophobic bilayer thickness (see Eq. 9), so it is not independent of . The second-order difference between and , which is shown in Eq. 10, arises from this difference in projection between actual thicknesses and equilibrium thicknesses within the definition of .)

On tension.

First of all, existing models [7], [12][15], [18], [22] were constructed at zero applied tension, which means in Eq. 2. To our knowledge, our work is the first where the coefficient of the term linear in is explicitly related to the applied tension (see Methods, Sec. 2) and to the tension of the Helfrich model (see Methods, Sec. 1.3).

In Ref. [18], the effect of applied tension is taken into account, in so far as it changes the equilibrium membrane thickness of a homogeneous membrane, but without being fully implemented in the elastic model. Our more complete description gives back this effect on membrane thickness, when it is applied to the particular case of a homogeneous membrane (see Methods, Sec. 2).

On the constant .

In our model, the constant features three contributions with different origins (see Eq. 4).

The first contribution arises from the spontaneous curvature of a monolayer and from its variation with the area per lipid. More precisely, the term(11)appears when one constructs the membrane model starting from a monolayer Hamiltonian density such as Eq. 1. This term was first introduced in Ref. [12], and it was then included in Refs. [13], [14].

The second contribution, , arises from , i.e., from the term in introduced in Eq. 1. This term was not included in Refs. [12][14], which started from a second-order expansion of the effective Hamiltonian per lipid molecule involving only the curvature and the area per lipid. However, a gradient of area per lipid (or, equivalently, of the thickness) in a monolayer has an energetic cost of its own. Indeed, a greater part of the hydrophobic chains is in contact with water when a thickness gradient is present (see Fig. 2). The associated energetic cost is given by the interfacial tension of the hydrocarbon-water interface, which is of order 40–50 mN/m [24], [25]. Such a term is often accounted for in microscopic membrane models (see, e.g., Ref. [26]). In the case of a symmetric membrane () with flat average shape, the surface of the hydrocarbon-water interface is increased by a factor for each monolayer (see Fig. 2). Thus, to second order, the associated energetic cost per unit projected area is . Note that other physical effects, e.g., the elasticity of the chains, may yield contributions to the term in . However, if we restrict to the simple term arising from interfacial tension, we obtain(12)

thumbnail
Figure 2. Thickness gradient.

Cut of a bilayer membrane with a symmetric thickness gradient. The dashed blue lines correspond to the hydrocarbon-water interfaces.

https://doi.org/10.1371/journal.pone.0048306.g002

Finally, the third contribution, , arises from the (macroscopic) externally applied tension. The tension of a vesicle can rise only up to a few mN/m before it bursts (see, e.g., Ref. [18]). Hence, according to our estimate of in Eq. 12, we expect .

In the seminal article Ref. [7], where the membrane model was constructed by analogy with liquid crystals, a term in , interpreted as arising from tension, was included in the effective Hamiltonian. However, its effect was neglected on the grounds that the value of its prefactor made it negligible with respect to the other terms. The value of this prefactor was taken to be that of the tension of a monolayer on the surface of a Plateau border [27]. The model introduced in Ref. [7] was further developed and analyzed in Refs. [18], [22], where the same argument was used to neglect the term in .

However, our construction of the membrane effective Hamiltonian shows that the microscopic tension involved through arises from local variations in the area per lipid. This stands in contrast with the case of the Plateau border, where whole molecules can move along the surface and exchange with the bulk, yielding a smaller value of the tension. Ref. [27] stresses that the measured tension of a Plateau border is valid for long-wavelength fluctuations, but that it is largely underestimated for short-wavelength fluctuations (less than 10 nm) which involve significant changes in area per molecule.

Including the tension of the hydrocarbon-water interface instead of that of the Plateau border is a significant change, given that the former is of order 40–50 mN/m [24], [25], while the latter is of order 1.5–3 mN/m [7], [18], [22], [27]. In Refs. [18], [22], it is shown that the effect of the term in is negligible if(13)where we have used our own notations of the prefactors of the terms in , and . In the case of DOPC, taking and using the values of the membrane constants [28], this condition becomes . While this is well verified if corresponds to the tension of the Plateau border, it is no longer valid within our model.

Our model is the first that includes all contributions to , in particular the one arising from interfacial tension. Besides, in Sec. 1.2 of our Methods part, we show that is also involved in , which emphasizes the complexity of constructing a continuum model to describe membrane elasticity at the nanoscale: many terms involved in the expansion of the effective Hamiltonian cannot be neglected a priori.

In the following, we will analyze numerical and experimental data, looking for evidence for the presence of , and comparing the relative weight of the different contributions to .

On the value of .

We have obtained (see Eq. 5), where is the bending rigidity of a symmetric membrane such that . The elastic constant is related to the bending rigidity of the Helfrich model (see Methods, Sec. 1.3) through(14)The difference between and arises from integrating out (see Methods, Sec. 1.2). In the previous models, this procedure was not carried out, as one focused directly on the symmetric case . All previous models thus made the approximation [7], [12][14], [18], [22].

In addition, in Sec. 1.2 of our Methods part, we show that is also involved in , which stresses further the possible importance of such terms in order to describe membrane elasticity at the nanoscale.

On boundary terms.

The boundary terms correspond to the last three terms in Eq. 2. When one wishes to describe the local membrane deformation due to a transmembrane protein, boundary terms play an important part, as their integral on the contour of the protein contributes to the deformation energy. The first two boundary terms are the same as in Refs. [12][14]. However, even at vanishing applied tension, we have , contrary to the previous models [14], due to the presence of . We have also accounted for the Gaussian bending rigidity , as in Ref. [15]: it yields the third boundary term.

Again, the situation is more complex when is included, as the expressions of and then feature extra terms linear in (see Eq. 37 in Sec. 1.2 of our Methods part).

On lipid tilt.

Several membrane models including lipid tilt in addition to average shape deformations and/or thickness deformations have been elaborated [21], [23], [26], [29][31]. These models provide improvements with respect to the Helfrich model, yielding better agreement with numerical data on bulk membranes [23], [31].

Our model does not include lipid tilt because we focus on local thickness deformations, and especially on comparison to experimental and numerical data regarding deformations induced by mismatched proteins. While it would be interesting to include this extra degree of freedom, it would imply introducing several membrane parameters, which would make comparison to mismatch data impractical.

Not taking tilt into account means that we are effectively integrating out this degree of freedom through coarse-graining. More precisely, the elastic coefficients of a more detailed membrane model, which would include tilt as an extra degree of freedom, would be renormalized by integrating out tilt. This means that tilt is included within the elastic coefficients of our membrane model. In addition, the interaction energy between the membrane and a mismatched inclusion (see, e.g., Eq. 15), and, consequently, the effective boundary conditions at the inclusion boundary, may involve tilt (see, e.g., Ref. [21]). In this interaction energy, tilt can be integrated out in the same way as in the bulk membrane energy. Hence, we are not losing any part of the elastic energy by disregarding the tilt degree of freedom. However, it is not impossible that a model including tilt truncated at second order could prove more efficient (e.g., have a wider domain of validity at short wavelengths) than one truncated at the same order and disregarding tilt.

Comparison with numerical results

As numerical simulations become more and more realistic, they start providing insight into the behavior of systems on the microscopic scale where direct experimental observation is difficult. Lipid membranes (with or without inclusions) are no exception. Over the last decade, several groups have simulated bilayer systems over length- and time-scales long enough to give access to the material constants relevant for nanoscale deformations. These simulations provide interesting tests for theoretical models describing membrane elasticity at the nanoscale. We will compare the predictions of our model to recent numerical results in this Section. All the numerical results we will discuss have been obtained at zero applied tension. Hence, throughout this section, we take . This implies that our definition of the membrane thickness is equivalent to that considered in the original numerical works (see the discussion above on the definition of ).

Fluctuation spectra.

Using numerical simulations, one can measure precisely the fluctuation spectra of the average height and the thickness of a bilayer membrane [14], [16], [32], [33]. Microscopic protrusion modes, occurring at the scale of a lipid molecule, contribute to these spectra. While they are not described by continuum theories, it is possible to consider that they are decoupled from the larger-scale modes [14], [16]. By fitting the numerical spectra to theoretical formulas, it is possible to extract the numerical values of the membrane constants involved in the continuum theory. In our framework, the fluctuation spectra of the average height of the membrane give access to the Helfrich bending rigidity , while those regarding the thickness of the membrane give access to , and .

We have reanalyzed the height and thickness spectra presented in Refs. [16], [32], [33] using the fitting formulas in Refs. [14], [16] (see Eq. 32 of Ref. [14]) and the method described in Ref. [14], except that we did not assume that , in order to include the possible effect of the difference between and (see Eq. 33), and of (see Eq. 35). Our results were similar to those obtained in Refs. [14], [16] assuming that , and we obtained no systematic significant difference between and , which means that the corrections to predicted by our model are negligible in these simulations. This gives a justification for focusing only on the correction to , as we do in this article. Besides, we obtained from all the fits, as reported in Refs. [14], [16], and we checked that all the values obtained for comply with the stability condition Eq. 8.

Deformation profiles close to a mismatched protein.

In Refs. [14][16], the thickness profile of a membrane containing one cylindrical inclusion with a hydrophobic mismatch has been obtained from coarse-grained numerical simulations. Comparing the average numerical thickness profiles to the equilibrium profiles predicted from theory is a good test for our model, in particular to find clues for the presence of .

Let us denote the radius of the protein by , and its hydrophobic length by : the mismatch originates from the difference between and the equilibrium hydrophobic thickness of the membrane. The equilibrium shape of the membrane, which minimizes its deformation energy, is solution to the Euler-Lagrange equation associated with the effective Hamiltonian density in Eq. 2. We write down this equilibrium shape explicitly in Sec. 3.1 of our Methods part. In order to determine it fully, it is necessary to impose boundary conditions at the edge of the inclusion, i.e., in . There is a consensus on the assumption of strong hydrophobic coupling , as it costs more energy to expose part of the hydrophobic chains to water than to deform the membrane, for typical mismatches of a few Å. Note that, with our definition of , the condition is valid to first order, while it is exactly valid with the definition of Refs. [7], [12][15], [18], [22], [23] (see Eqs. 9, 10). This difference arises from the fact that our is not projected along (see Fig. 1), which makes it fully independent of . Given that the elastic energy is known to second order, the equilibrium membrane shape resulting from its minimization is known to first order, so it is sufficient to use boundary conditions to first order. Hence, such differences are not relevant for the present study and will not be mentioned any longer.

However, there is some debate about the second boundary condition in (see, e.g., Ref. [14]), which regards the slope of the membrane thickness profile. Traditionally, one either assumes that the protein locally imposes a fixed slope to the membrane [18], [22], or minimizes the effective Hamiltonian in the absence of any additional constraint, which amounts to considering that the system is free to adjust its slope in [12][16]. In Sec. 3.1 of our Methods part, we present the equilibrium profiles for these two types of boundary conditions. The actual boundary condition depends on the interactions between the protein and the membrane. In a quadratic approximation, these interactions generically give rise to an effective potential favoring a slope in :(15)where is an effective rigidity, while denotes the derivative of the membrane thickness profile with respect to the radial coordinate . Two a priori unknown parameters, and , are associated with this effective potential. The “free-slope” boundary condition (also called “natural” boundary condition [12], [14]) is recovered in the limit , which is appropriate if is negligible with respect to the energetic contributions in . Conversely, if , the protein locally imposes the fixed slope . If the interactions between the protein and the membrane lipids are sufficiently short-ranged, the protein cannot effectively impose or favor a slope on the coarse-grained membrane thickness profile. For instance, in the numerical simulations of Refs. [14][16], the interactions between the protein and the membrane lipids are of similar nature and of similar range as those between membrane lipids. Thus, we will choose the free-slope boundary condition in our analysis of this data. This choice was already made in Refs. [14][16]. A practical advantage of this boundary condition is that it does not introduce any unknown parameter in the description.

The membrane model of Refs. [14][16] is very similar to ours, except that . It was shown in Ref. [16] that this model can reproduce very well the numerical results, provided that the spontaneous curvature is adjusted for each deformation profile (see Fig. 3). In Ref. [16], the adjusted “renormalized spontaneous curvature”, denoted by , was found to depend linearly on the hydrophobic mismatch [16], as shown in Fig. 4. In our model, the equilibrium profile corresponding to the free-slope boundary conditions (see Eqs. 46 and 53) involves . We show in Sec. 3.1 of our Methods part that the quantity(16)then plays the part of the renormalized spontaneous curvature of Ref. [16] in the equilibrium profile. This quantity is linear in : our model, and more precisely the presence of a nonvanishing , thus provides an appealing explanation for the linear dependence observed in Ref. [16].

thumbnail
Figure 3. Thickness deformation due to a mismatched inclusion.

Membrane thickness profile from Ref. [16] in the vicinity of a mismatched inclusion with hydrophobic thickness and radius , with center in , as a function of the radial coordinate . The equilibrium membrane hydrophobic thickness is . The unit of length on the graph is 6 Å, as in Ref. [16]. Dots: numerical data (the error bars on the data, not reproduced here, are about 1 Å wide [16]). Line: best fit. Exactly as in the original reference, the numerical data is fitted to Eqs. 4653 with , taking and the (renormalized) spontaneous curvature as fitting parameters, the other constants being known from the fluctuation spectra.

https://doi.org/10.1371/journal.pone.0048306.g003

thumbnail
Figure 4. Renormalized spontaneous curvature

as a function of the hydrophobic mismatch . Data from Ref. [16], which presents fits of simulation results for inclusions with three different hydrophobic thicknesses. Line: linear fit, with slope . Note that our corresponds to twice that in Table 2 of Ref. [16], as we work with total curvatures instead of average curvatures. The error bars on are those listed in that table, and corresponds to in that table.

https://doi.org/10.1371/journal.pone.0048306.g004

Using a linear fit of the data of Ref. [16] (see Fig. 4), together with Eq. 16 and the value extracted from the spectra in Ref. [16], we obtain .

It is interesting to compare this value to , which is obtained from the fluctuation spectra in Ref. [16]: . This shows that the contribution of to is important. Besides, we may now estimate the contribution to that stems from the monolayer spontaneous curvature (see Eq. 4): . Using values from the fluctuation spectra in Ref. [16], this yields for the algebraic distance from the neutral surface of a monolayer to the hydrophilic-hydrophobic interface of this monolayer (see Methods, Sec. 4 for the relation between and ).

In Ref. [15], a different coarse-grained molecular simulation model was used to obtain the equilibrium membrane thickness profiles for cylindrical inclusions with two different hydrophobic thicknesses, one yielding a positive mismatch and the other a negative one, and with seven different radii . These profiles are presented in Figs. 6 and 7 of Ref. [15], except those corresponding to the inclusions with largest radii (5.25 nm), but this data was communicated to us by one of the authors of Ref. [15]. We fitted each of these numerical profiles to the analytical equilibrium profile Eq. 46 with prefactors (see Eq. 54), using as our only fitting parameter, in the spirit of Ref. [16]. We found that does not depend on the radius of the inclusion, but that it depends significantly on the mismatch (see Fig. 4A). This is in good agreement with the predictions of our model (see Eq. 16). For each of the two values of , we have averaged over the seven results corresponding to the different inclusion radii. The line joining these two average values of as a function of is plotted in Fig. 5B. Using Eq. 16 and the value [14], [15], the slope of this line yields : the order of magnitude of this value is the same as the one obtained from the data of Ref. [16].

thumbnail
Figure 5. Renormalized spontaneous curvature

as a function of the inclusion radius and the hydrophobic mismatch . A) versus . The values of were obtained by fitting each thickness deformation profile of Ref. [15]. Circles (blue): positive mismatch, . Squares (red): negative mismatch, . Solid lines: average values; dotted lines: standard deviation over the seven data points (corresponding to the different ), for each value of . B) Average value of (see A) as a function of the hydrophobic mismatch . The equation of the line joining the two data points has a slope .

https://doi.org/10.1371/journal.pone.0048306.g005

Again, we can compare this value to , which is obtained from the fluctuation spectra in Refs. [14], [15]: . Hence, the contribution of to is important here too. We also obtain , and .

In Ref. [15], the shortcomings of the model that disregards are explained by the local variation of the volume per lipid close to the protein. It is shown in Ref. [15] that including this effect yields(17)where is the bulk equilibrium volume per lipid, while denotes the volume per lipid in . However, the predicted linear dependence of in is not obvious: in Fig. 6, we rather see two groups of points (one for each value of ) than a linear law. In other words, the data of Ref. [15] is more consistent with a value of that depends only on and not on (or ), in agreement with the predictions of our model (see Eq. 16). In Ref. [16], local modifications of the volume per lipid close to the inclusion were investigated too, as well as local modifications of the nematic order, of the shielding of the hydrophobic membrane interior from the solvent, and of the overlap between the two monolayers. None of these effects was found to explain satisfactorily the linear dependence of versus [16].

thumbnail
Figure 6. Renormalized spontaneous curvature

versus the relative volume variation on the inclusion edge. The values of are extracted from fitting the data of Ref. [15], and the values of are directly taken from Ref. [15].

https://doi.org/10.1371/journal.pone.0048306.g006

To sum up, our model can explain the dependence of in observed in the numerical results of Refs. [15], [16] as a consequence of the presence of . Our explanation does not involve any local modification of the membrane properties, in contrast with those proposed in Refs. [15], [16]. Furthermore, the order of magnitude we obtain for from the data of Refs. [15], [16] is in agreement with our estimate in Eq. 12.

Comparison with experimental results

The antimicrobial linear pentadecapeptide gramicidin (see [8] for a review) is a very convenient experimental system to probe membrane elasticity on the nanoscale. In lipid membranes, two gramicidin monomers (one in each monolayer) associate via the N-terminus to form a dimeric channel, stabilized by six intermolecular hydrogen bonds. The channel being large enough for the passage of monovalent cations, conductivity measurements [9] can detect its formation and lifetime, which are directly influenced by the membrane properties. Indeed, while the monomers do not deform the membrane, the dimeric channel presents a hydrophobic mismatch with the membrane, so that dimer formation involves a local deformation of the bilayer. The gramicidin channel can therefore act as a local probe for the bilayer elasticity. Furthermore, the gramicidin channel can be considered as up-down symmetric and cylinder-shaped, which makes it convenient for theoretical investigations.

Data on gramicidin channels originally motivated theoretical investigations on membrane models describing local thickness deformations [7], [11][13]. Such data now provides a great opportunity to test any refinement of these models. We will compare our model to the data of Ref. [17] regarding the lifetime of the gramicidin channel as a function of bilayer thickness, and then to the data of Ref. [18] regarding the formation rate of the gramicidin channel as a function of bilayer tension.

In order to compare the predictions of our model to experimental data regarding the gramicidin channel, it is necessary to make some assumptions about the boundary conditions at the edge of the channel, i.e., in . As discussed in the previous section, we will assume strong hydrophobic coupling, i.e., , but determining the boundary condition on the slope of the membrane thickness profile is trickier as it depends on the interactions between gramicidin and the membrane lipids. In previous analyses [18], [34], the fixed-slope boundary condition was favored as giving the best agreement with experimental data. However, different values of the fixed slope were obtained in these studies. In addition, recent all-atom simulations of gramicidin channels in lipid bilayers indicate that the membrane thickness profile is complex in the first lipid shell around the channel, due to specific interactions, and that beyond this first shell, no common slope exists for the different membranes investigated [35]. Given the difficulty to determine the actual effective boundary condition associated with the slope of the membrane thickness profile, we will adopt the free-slope boundary condition, which has the advantage not to introduce any unknown parameter in the analysis, but we will also compare our results to those obtained with the more traditional fixed-slope boundary condition.

Analysis of the experimental data of Elliott et al. [17].

It was shown in Ref. [22] that the detailed elastic membrane model introduced in Ref. [7] yields an effective linear spring model as far as the membrane deformation due to gramicidin is concerned [22], [34]: the energy variation associated with the deformation can be expressed as , where is the effective spring constant, while is the thickness mismatch between the gramicidin channel and the membrane. This linear spring model was validated by comparison with experimental data on the lifetime of the gramicidin channel, measured as a function of bilayer thickness ([17], [36], summarized in [34]) and as a function of the channel length [37].

We will here focus on the data concerning virtually solvent-free bilayers, i.e., membranes formed using squalene. The elasticity of membranes containing hydrocarbons should be different: for instance, a local thickness change of the membrane could be associated with a redistribution of the hydrocarbons. (In this, our analysis differs from that of Ref. [14], where all the data of Ref. [17] was considered. Another important difference with the analysis conducted in that reference is that we use experimental values of the membrane parameters, which are quite different from the values coming from numerical simulations.) In Ref. [34], the effective spring constant of the membrane was estimated from data of Ref. [17] on gramicidin channel lifetime for three bilayers formed in squalene with monoglycerids that differed only through their chain lengths: the different thicknesses of these membranes yield different hydrophobic mismatches with a given type of gramicidin channels. The value was obtained.

In Sec. 3.2 of our Methods part, we use our model to calculate the deformation energy of the membrane due to the presence of a mismatched protein. Both in the case of the free-slope boundary condition, and in the case where the gramicidin channel locally imposes a vanishing slope, this deformation energy can be expressed as a quadratic function of the mismatch . The prefactor of in the deformation energy corresponds to the effective spring constant of the system. Thus, although our model is different from the one of Refs. [7], [18], [22], it also yields an effective linear spring model. This is not surprising since we are dealing with the small deformations of an elastic system. However, the detailed expressions of our spring constants as a function of the membrane parameters (see Eqs. 59 and 65) are different from those obtained using the model of Refs. [7], [18], [22], due to the differences between the underlying membrane models. In particular, in our model, is involved in , through . Our aim will be to find out which value of gives the best agreement with the experimental value of .

Using Eqs. 4, 5 and 7, and neglecting the difference between and , Eqs. 59 and 65 show that depends on the elastic constants , and involved in the Helfrich model, on , on , which corresponds to the spontaneous curvature variation with the area per lipid, on , on the radius of the gramicidin channel, and on . There is, to our knowledge, no direct experimental measurement of available, but, as shown in Sec. 4 our Methods part, we have , where denotes the algebraic distance from the neutral surface of a monolayer to the hydrophilic-hydrophobic interface of this monolayer (see Eq. 73, neglecting the difference between and ). Hence, in order to calculate the spring constant, we need values for , , , and , in the precise case of monoolein membranes.

In Ref. [38], the elastic constants , and were measured in a monoolein cubic mesophase, both at and at . The positions of the neutral surface and of the hydrophilic-hydrophobic interface were estimated on the same system in Ref. [39], but these results were flawed by a mathematical issue, which was corrected in Ref. [40]. This correction yielded other corrections on , and on the ratio [41]. These results regard a cubic phase, where the membrane is highly deformed with respect to a flat bilayer: the values of the various constants should be affected by the strains present in this phase. In another work [42], the constants of monoolein are determined in a highly hydrated doped phase, where the strains are better relaxed. However, these measurements were carried out at , while the experiments of Ref. [17] that we wish to analyze were performed at . Given that the data of Refs. [38], [39] include the most appropriate temperature, while the ones of Ref. [42] correspond to the most appropriate phase, we will present results corresponding to both sets of parameters. Finally, the experimental value of for monoolein is provided by Ref. [27].

In Table 1, we present the results obtained for the spring constant of monoolein bilayers, using the different experimental estimates of the membrane constants. The main difference between parameter sets 1 and 2 is the value and the sign of [38], [41]. However, is involved in only in the free-slope case (see Eqs. 59 and 65): the 3% difference between the values of obtained with parameter sets 1 and 2 stems only from the difference on , while the 12% difference between obtained with data sets 1 and 2 contains an important contribution from . The constants in parameter set 3, corresponding to Ref. [42], are significantly different from those of Refs. [38], [41], which yields a 30% difference on and a 20% difference on . We also note that, as the value of the algebraic distance from the neutral surface to the hydrophilic-hydrophobic interface of a monolayer is very small compared to the other length scales involved ( [40]), the contribution of to is negligible (it is of order 1%).

Let us now discuss the results given by our model, in the case of the free-slope boundary condition (see Table 1). The spring constants obtained assuming that are about three times smaller than the experimental value (see line 1 of Table 1). (This result is very similar to that in Ref. [34], which illustrates that accounting for monolayer spontaneous curvature and for boundary terms does not change much the value of .) However, reaches the experimental value for for all three parameter sets (see line 2 of Table 1). Hence, for free-slope boundary conditions, the presence of , with an order of magnitude consistent with Eq. 12, improves the agreement between theory and experiment.

We may compare these values of to the contribution to that originates from the monolayer spontaneous curvature (see Eq. 4): . We estimate the value of this contribution to be between and , depending on which set of parameters is chosen. This is positive and much smaller in absolute value than the estimates obtained from the numerical data of Ref. [16] and of Ref. [15]: here, the neutral surface of a monolayer and its hydrophilic-hydrophobic interface are very close, while seemed to be significant in the numerical simulations. In addition, the contribution of membrane tension to , namely, , cannot exceed about 1 mN/m. In the case of the free-slope boundary condition, our results imply that should be the dominant contribution to for the membranes studied in Ref. [17].

Let us now discuss the results obtained for the zero-slope boundary condition, which was investigated in Ref. [34]. For the zero-slope boundary condition, the values obtained for assuming that are in quite good agreement with the experimental value obtained in Ref. [34] from the data of Ref. [17], for all the data sets we used (see line 3 of Table 1): hence, seems negligible if zero-slope boundary conditions are assumed. However, there is no justification to assume that the gramicidin channel locally imposes a vanishing slope.

Analysis of the experimental data of Goulian et al. [18].

While the experiments cited in the previous Section dealt with discrete changes of the hydrophobic mismatch obtained by varying membrane composition, Goulian et al. [18] measured the gramicidin channel formation rate in lipid vesicles as a function of the tension applied through a micropipette. As the tension is an externally controlled parameter that can be changed continuously for the same gramicidin-containing membrane, this approach can yield more information, and it has the advantage of limiting the experimental artifacts associated to new preparations. To date, the experiment in Ref. [18] remains the most significant in the field and should serve as a testing ground for any theoretical model. We will therefore discuss in detail the data and its interpretation by the original authors [18], [22] as well as in terms of our model (see Eq. 2).

Within experimental precision, the data of Ref. [18] can be described by a quadratic dependence:(18)Given that is a linear function of the energy barrier associated with the formation of the gramicidin dimer, it is a sum of a chemical contribution, including, e.g., the energy involved in hydrogen bond formation, and of a contribution arising from membrane deformation due to the dimer (monomers do not deform the membrane) [18]. The latter contribution arises from the hydrophobic mismatch between the membrane and the dimer, and it depends on the applied tension , since the membrane hydrophobic thickness depends on (see Eq. 43 in Sec. 2 of our Methods part). Expressing the deformation energy of the membrane due to the presence of the dimer gives a theoretical expression for the -dependent part of . In our model, features a contribution coming from (see Eq. 4). However, this term is negligible, given that (see Eq. 13), for realistic tension values (a few mN/m at most), and for the experimentally measured values of the membrane constants [28]. This enables us to disregard it. Then, our quadratic elastic membrane model simply gives a quadratic dependence of on , in agreement with the form of Eq. 18. Comparing the experimental values of and to those predicted by theory provides a test for theoretical models [18]. (Note that, if the -dependent contribution to is included, the expression of the -dependent part of is no longer simply quadratic in . However, we explicitly verified that including this contribution yields a negligible change to the relation between and , for realistic values of the parameters).

Since the coefficients and arise from membrane elasticity, they are common to all the vesicles studied in Ref. [18], which have the same lipid composition. Conversely, the baseline depends on parameters such as the concentration of gramicidin molecules, so it can take a different value for each of the twelve vesicles studied in Ref. [18]. A global fit to the data of Ref. [18] using Eq. 18 involves minimizing the goodness-of-fit function(19)where the index runs over all the experimental points, with fitting parameters . The baseline is then subtracted from each of the twelve curves. All the data is plotted in the same graph in Fig. 7. The best global fit, corresponding to and , is shown on Fig. 7 as the dotted (black) line. (It should be noted that the values obtained by fitting the individual curves are much more scattered: ranges from to and from to .)

thumbnail
Figure 7. Formation rate

of gramicidin channels versus the applied tension , analyzed with a quadratic model. Diamonds: experimental data retrieved from Fig. 6b of Ref. [18], after subtraction of the baselines . Dotted black line: best quadratic fit, with and ; . Dashed green line: results obtained from the elastic model of Ref. [22], with the constants given in [18]; . Dashed-dotted blue line: idem with more recent values of the constants; . Solid red line: results obtained by taking and the recent values of the constants in the model of Refs. [7], [22]; . The values of and corresponding to the curves on this graph are listed in Table 2.

https://doi.org/10.1371/journal.pone.0048306.g007

In Ref. [18], the authors used published values of the material constants to calculate and in the framework of their elastic model [22], based on that of Ref. [7]. Using fixed-slope boundary conditions, they reported good agreement with the experimental data for a reasonable value of the unknown slope (). However, we need to raise the following points:

  1. There was a mistake in their implementation of the formula of Ref. [22] giving and as a function of the material constants. More precisely, we found that a factor of 2 was missing in the expression of and a factor of 4 was missing in that of in the implementation of the formula of Ref. [22]. This was confirmed by Mark Goulian (private communication). The actual values of and obtained using the same values of the constants as in Ref. [18] are in fact quite far from those corresponding to the best fit of the experimental data, as shown by the dashed green line in Fig. 7 (see also Fig. 8 and Table 2).
  2. The estimates for the elastic constants used in Ref. [18] are somewhat different from more recent and more widely accepted values. Henceforth, we will use the following parameters, for a DOPC membrane: [18], , [28], [43], and the dimensions of a gramicidin channel: , [18]. Implementing these more recent values in the model of Ref. [22] does not yield a better agreement with experiment, as shown by the dashed-dotted (blue) line in Figure 7 (see also Fig. 8 and Table 2).
thumbnail
Figure 8. Comparison between the experimental values of

and and those obtained from different models. Colorscale: goodness-of-fit function (see Eq. 19) for the data of Ref. [18], as a function of the fitting parameters and . White diamond: values of and that give the best fit. Black triangle: results obtained from the elastic model of Ref. [22], with the constants given in [18]. Lines: trajectories obtained from our model in the plane when varying . Red: free slope; green: , black: . These three curves start by a white dot at , and increases rightwards along these curves. The rightmost white dot (, ) roughly corresponds to the best agreement we can obtain between our model and the experiment fitted to the quadratic model (red curve on Fig. 7). The black diamond corresponds to the best agreement we can obtain between our model and the experiment fitted to the linear model at low tension (see Fig. 9).

https://doi.org/10.1371/journal.pone.0048306.g008

thumbnail
Table 2. Values of and obtained from the model of Ref. [22] and from our model with .

https://doi.org/10.1371/journal.pone.0048306.t002

A somewhat better agreement with the experimental data is obtained when taking instead of for the fixed slope (see Figs. 7 and 8, and Table 2). However, the downward inflection of the experimental curves at high is not adequately described for any value of . In fact, is independent of , and its absolute value given by the elastic model is 15 times smaller than the experimental one (see Table 2). We conclude that the elastic model of Refs. [7], [22] does not satisfactorily describe the data of Ref. [18] regarding the lifetime of the gramicidin channel under tension.

In Sec. 3.2 of our Methods part, we calculate the deformation energy in the framework of our model, both for the fixed-slope boundary condition and for the free-slope boundary condition. The resulting expressions of and are given by Eqs. 61, 62, 67 and 68. In order to see which values of and which boundary conditions give the best agreement with the experiments of Ref. [18], we present a plot of the goodness-of-fit function (see Eq. 19) in a graph in Fig. 8. On this graph, we have plotted the trajectories obtained from our model in the plane when varying , for , for (as in Ref. [18]), and for the free-slope boundary condition.

In order to obtain numerical values of and from Eqs. 61, 62, 67 and 68, we used the above-mentioned parameter values, and the estimate [16]. Finally, we estimated through the relation (see Eq. 73 in Sec. 4 of our Methods part). For this, the algebraic distance from the neutral surface of a monolayer to the hydrophilic-hydrophobic interface of this monolayer was estimated by first determining the position of the pivot surface from the data of Ref. [43], and by calculating the distance between it and the neutral surface [44]: we found . Here again, the neutral surface is close to the hydrophilic-hydrophobic interface. For the sake of simplicity, we took , and we checked that the results were not significantly different when taking .

The ingredient in our model that can change significantly the results is (Note that the values of and corresponding to are very close to those obtained using the model of Ref. [18] with our values of the parameters, as shown in Table 2. This illustrates again that the influence of boundary terms is quantitatively small.) Fig. 8 shows that the experimental value of can be explained by our model. In addition, the values of that minimize , i.e., that give the best agreement with the experimental data of Ref. [18], are between and , depending on the boundary condition chosen, as shown in Table 3. This range of values of is reasonable.

thumbnail
Table 3. Values of , and obtained from our model that yield the best agreement with the experimental results of Ref. [18], analyzed with a quadratic fit (see Eq. 18 and Fig. 7).

https://doi.org/10.1371/journal.pone.0048306.t003

For the free-slope boundary condition, the best agreement with the experimental results is obtained for (see Table 3 and Fig. 8). The order of magnitude is the one expected from .

Let us now discuss the results obtained for the fixed-slope boundary condition, which is used in Ref. [18]. For a fixed slope , the best agreement with the results of Ref. [18] analyzed with the complete quadratic fit is obtained for . Conversely, for , the best agreement is obtained for , which is similar to the result obtained the free-slope case (see Table 3 and Fig. 8). Hence, in the case of the fixed-slope boundary condition, the conclusions depend a lot on the value of that is chosen.

In all cases, the absolute values of we obtain remain much smaller than the one that matches best the experimental results, which is (see Fig. 7). This can be seen in Fig. 8, as well as in Table 3. Hence, with our model, as with the one of Ref. [18], it seems impossible to explain the experimental value of . Our model predicts that is proportional to the effective spring constant of the membrane discussed in the previous Section (see Eqs. 59 and 65): it is thus quite unexpected to have a good agreement with the experimental values of but not with those of . This disagreement on could come either from a shortcoming of the model or from an undetected systematic error in the experimental data. The importance of is largest at highest tensions, as it is which gives the curve its concavity, and it should be noted that the maximum applied tension is around in Ref. [18], which is comparable to the rupture threshold of [28]. The membrane properties may be affected at such high tensions in a way that is no longer well described by standard elastic models. It would be interesting to have more experimental data on the behavior of gramicidin channels under tension to see if this unexpected value of persists.

Following the hypothesis that high tensions are problematic, we performed a linear fit of the data of Ref. [18] (i.e., a fit with ), keeping only the points corresponding to : this yields (see Fig. 9). In Table 4, we list, for different boundary conditions, the value of which gives , and the value of obtained from our model for this . These values correspond to those that give the best agreement between our model and the linear fit to the low-tension data of Ref. [18] presented in Fig. 9. Table 4 shows that the values of that yield the best agreement with the experimental data have a similar order of magnitude as those obtained above with the full quadratic fit (see Table 3), remaining below . Again, these values depend a lot on for fixed-slope boundary-conditions. (For instance, the slope is consistent with (see Table 4). However, there is no a priori reason for assuming that .)

thumbnail
Figure 9. Formation rate

of gramicidin channels as a function of the applied tension , analyzed with a linear model, for . Diamonds: experimental data retrieved from Fig. 6b of Ref. [18], after subtraction of the baselines (which are different from those of Fig. 7 since the fitting model is here linear instead of quadratic). Line: best linear fit, yielding ; correlation coefficient: .

https://doi.org/10.1371/journal.pone.0048306.g009

thumbnail
Table 4. Values of and obtained from our model that yield the best agreement with the experimental results of Ref. [18] analyzed with the low-tension linear fit.

https://doi.org/10.1371/journal.pone.0048306.t004

Again, we may compare our estimates of (see Tables 3 and 4) to the term , which also contributes to : here, . This is much smaller in absolute value than the corresponding estimates obtained from the numerical data of Ref. [16] and of Ref. [15]: here, as in the membranes studied in Ref. [17], the neutral surface of a monolayer and its hydrophilic-hydrophobic interface are very close, while seemed to be of a few Å in the numerical simulations. We note in passing that this hints at a relevant difference between simulated membranes and real membranes. Besides, in the case of the free-slope boundary condition, our results imply that should be the dominant contribution to for the membranes studied in Ref. [18], as for those of Ref. [17].

Hence, for the free-slope boundary condition, our analyses of the numerical data of Ref. [16] and of Ref. [15], and our analyses of the experimental data of Ref. [17] and of Ref. [18] all converge toward a value of a few tens of mN/m for , which is of the order of magnitude expected if . Conversely, for the fixed-slope boundary condition, the value of is coupled to that of the slope .

Conclusion

We have put forward a modification of membrane elastic models used to describe thickness deformations at the nanoscale. We have shown that terms involving the gradient (and the Laplacian) of the area per lipid contribute to important terms of the effective Hamiltonian of the bilayer membrane. We have reanalyzed numerical and experimental data to find some signature of the presence of these terms. Using the free-slope boundary condition at the boundary of the mismatched protein, we have obtained consistent results showing that the term stemming from the gradient of the area per molecule has a prefactor in the range . Such values are consistent with the idea that this term involves a significant contribution of the interfacial tension between water and the hydrocarbon-like hydrophobic part of the membrane. Indeed, this contribution should yield .

Interestingly, our analysis of the experimental data from Ref. [18] has shown that these nice experimental results were not as well understood as assumed in the literature. Hence, it would be interesting to have more data on the behavior of gramicidin channels in membranes under tension.

Finally, the effective linear spring model [22], [34] is a very useful simplification of membrane elastic models when dealing with local thickness deformations and hydrophobic mismatch. Its applicability has been thoroughly tested on systems where gramicidin is used to probe the influence of various molecules on membrane properties (see, e.g., Ref. [10]). As other quadratic elastic models, our model yields an effective spring model. However, since the expression of the spring constant depends on the details of the model, careful consideration is required when one is interested in the behavior of a particular material constant.

Methods

1 Derivation of the effective Hamiltonian

1.1 General expression of the bilayer effective Hamiltonian.

Let us consider a patch of bilayer membrane with a fixed projected area , at fixed chemical potential . The rest of the membrane (e.g., of the vesicle) plays the part of the reservoir that sets the chemical potential . The effective Hamiltonian per unit projected area in each monolayer is , where is given by Eq. 1, while the projected area per molecule reads to second order. Hence, Eq. 1 yields, to second order in the deformation and in the relative stretching of the monolayers,(20)

We assume that the hydrophobic chains of the lipids are incompressible. Let us introduce the excess hydrophobic thickness (resp. ) of the upper (resp. lower) monolayer, defined as its hydrophobic thickness along the normal to its hydrophilic-hydrophobic interface minus the equilibrium monolayer hydrophobic thickness (see Fig. 1). In the spirit of Refs. [12][14], we use the incompressibility condition(21)where is the constant hydrophobic volume per lipid. (In this incompressibility condition, a correction arising from membrane curvature is neglected. Using the complete incompressibility condition instead of this one yields the same effective Hamiltonian Eq. 2, but with different expressions of and as a function of the constants involved in Eq. 1. These expressions depend on , and consequently on the applied tension, but this dependence is negligible for realistic tension values. As the rest of our discussion is not affected by this, we keep the approximate incompressibility condition for the sake of simplicity. Note that the exact incompressibility condition was implemented recently in Ref. [23].)

In all the following, we will work to second order in the small dimensionless variables , , , and . In this framework, using the relations and , Eq. 20 becomes(22)In this expression, we have introduced the constitutive constants of a monolayer: is compressibility modulus of the monolayer, is its bending rigidity, is its Gaussian bending rigidity, is its spontaneous (total) curvature, and is the modification of the spontaneous (total) curvature due to area variations. More precisely, where is the lipid area-dependent (total) spontaneous curvature of the monolayer. In addition, recall that denotes the equilibrium hydrophobic thickness of the bilayer membrane. Finally, we have introduced the constants(23)(24)These two constants have the dimension of a surface tension, like .

In our description, the state of monolayer is determined by the two variables and . Hence, the state of the bilayer membrane is a priori determined by four variables. However, given that there must be no space between the two monolayers, the distance along between the hydrophilic-hydrophobic interfaces of the two monolayers must be equal to the sum of their projected thicknesses. Hence, to second order, we have the following geometrical constraint:(25)This leaves us with only three independent variables to describe the state of the membrane. Let us choose the average shape of the bilayer, the sum of the excess hydrophobic thicknesses of the two monolayers, and the difference between them:(26)(27)(28)

Thus, we can rewrite the effective Hamiltonian per unit projected area of the membrane in terms of the new variables , and . It reads, to second order in the small dimensionless variables , , , , , , , and , and discarding derivatives of order higher than two:(29)where we have introduced , which plays the part of an externally applied tension (see Methods, Sec. 2).

1.2 Eliminating .

In the present study, we are not interested in the variable . In a coarse-graining procedure, this degree of freedom can be eliminated by integrating over it. In our Gaussian theory, it simply amounts to minimizing with respect to . This variable is coupled to the membrane curvature , but not to . In the case of a constant curvature, the constant value(30)is a simple solution to the Euler-Lagrange equations in , for which the term involving in reads(31)As the variable varies spontaneously on length scales much shorter than the variable , we can consider in a first approximation that will simply follow , in which case this constant solution is the valid one. Thus, after this partial minimization, this term provides a correction to .

We finally obtain(32)where the usual Helfrich bending rigidity , associated with the average shape, is related to through(33)

In the case where the average shape of the membrane is flat, i.e., , dropping constant terms, we obtain the expression of in Eq. 2 with(34)(35)(36)(37)

Thus, in general, in Eq. 2, the constants , include contributions in and , which arise from , and (see Eqs. 23, 24). Therefore, the terms in gradient and Laplacian of introduced in Eq. 1 cannot be neglected a priori, as they contribute to the terms in and that are traditionally accounted for in models describing membrane thickness deformations [7], [12][14], [18], [22]. Due to these contributions, the values of the constants and are not fully predicted by the constants involved in the Helfrich model. This stands in contrast with the models developed previously [7], [12][14], [18], [22]. In addition, the terms arising from , and modify the relations between the various coefficients: in the previous models that accounted for boundary terms, assuming , and disregarding tension, one had [14], which is no longer true here. This will affect the equilibrium thickness profile of a membrane containing a mismatched protein.

1.3 Link with the Helfrich Hamiltonian.

Since the variables and are decoupled in the Hamiltonian density given by Eq. 32, the terms depending on can be isolated, yielding(38)which corresponds to the Helfrich Hamiltonian [3] for a membrane composed of two identical monolayers. In particular, the term in has the standard form of a Helfrich tension term, conjugate to the actual area of the membrane, since the element of area is to second order. Hence, can be viewed as an effective applied tension. This interpretation of is explained in more detail in Sec. 2 of our Methods part.

Hence, our model gives back the Helfrich Hamiltonian if the state of the membrane is described only by its average shape , i.e., if the variable is integrated out.

1.4 Stability criterion.

Let us focus on a membrane with flat average shape , described by Eq. 2. Depending on the values of the constants , and , a homogeneous thickness can be less or more energetically favorable than an undulated shape. The physical situation we wish to describe is the one where the equilibrium state has a homogeneous thickness. To determine which sets of constants comply with this, let us calculate the effective Hamiltonian per unit projected area of a membrane with harmonic undulations characterized by the wave vector . Neglecting boundary terms (by taking appropriate boundary conditions or by assuming that the undulations decay on some large length scale), we obtain , where the omitted prefactor is positive. The flat shape is favored if for all , and otherwise there exist some values of for which it is unstable. Thus, the conditions for the stability of the flat shape are , and .

2 Membrane submitted to an external tension

In Sec. 1 of our Methods part, we have derived the effective Hamiltonian of a bilayer membrane in the ensemble. This is the most convenient thermodynamic ensemble to work in. However, in order to describe experiments where a vesicle is submitted to an external tension, one should work in the ensemble, where is the number of lipids in the vesicle and is the externally applied tension. This is especially interesting in order to analyze the results of Ref. [18]. The ensemble change can be performed using a Legendre transformation: in the ensemble, the adapted effective Hamiltonian is , where , with expressed in Eq. 32, and(39)

Let us restrict ourselves to the case of a homogeneous and flat membrane, i.e., to a membrane with constant and . Then, using Eq. 39 to eliminate the variables and from the expression of , we obtain, to second order:(40)

Minimizing with respect to yields the equilibrium excess thickness of the membrane at a given imposed tension . To first order, it reads(41)Note that, since is assumed to be a first-order quantity, must be first-order too for our description to be valid for . This property has been used to simplify the result in Eq. 41. In practice, is well verified, given that cannot exceed a few mN/m without the vesicle bursting, while is of order . Since is the equilibrium hydrophobic thickness of this piece of homogeneous and flat membrane submitted to a vanishing external tension, it is consistent that vanishes when does, as is the excess thickness with respect to . Eq. 41 shows that the thickness of a membrane with fixed number of lipids decreases when the external tension increases, and is in agreement with Ref. [18].

We are now going to show that the constant in the ensemble (see, e.g., Eq. 32) plays the part of an externally applied tension. For this, let us calculate the equilibrium thickness of a membrane patch with projected area at a chemical potential , when it is homogeneous and flat. This amounts to minimizing with respect to . For a homogeneous and flat membrane, Eq. 32 becomes(42)Minimizing with respect to then gives(43)Comparing Eq. 43 to Eq. 41 shows that plays the part of the externally applied tension . Hence, can be considered as an effective applied tension.

3 Membrane containing a cylindrical mismatched protein

In this Section, we write down explicitly the equilibrium shape and the deformation energy of a membrane which contains a single cylindrical transmembrane protein with a hydrophobic mismatch (see Fig. 1B). This protein can correspond to a gramicidin channel in the dimer state. We focus on a membrane with a flat average shape, described by the effective Hamiltonian per unit projected area in Eq. 2. We denote the radius of the protein by , and its hydrophobic thickness by . We take the center of the cylindrical protein as the origin of the frame, which yields cylindrical symmetry.

In order to treat the case where the membrane is submitted to a tension , we rewrite Eq. 2 in terms of the variable , which represents the excess hydrophobic thickness of the bilayer relative to its equilibrium value at an applied tension (see Eq. 43). Discarding constant terms and using the relation , which yields , it yields(44)

3.1 Equilibrium thickness profile.

Let us first review (see, e.g., Ref. [22]) the equilibrium thickness profile of the membrane containing the mismatched protein. This equilibrium shape is solution to the Euler-Lagrange equation associated with the effective Hamiltonian in Eq. 44,(45)

Using the cylindrical symmetry of the problem and choosing solutions that vanish at infinity, we obtain, if the stability condition Eq. 8 is verified, the following solution to the Euler-Lagrange equation Eq. 45:(46)where is the -order modified Bessel function of the second kind, and(47)which are either both real or complex conjugate.

The integration constants are determined by the boundary conditions at . The first boundary condition corresponds to strong hydrophobic coupling: on the inclusion boundary, the hydrophobic thickness of the membrane is equal to that of the inclusion, which is denoted by (see Fig. 1B). It yields (to first order, as explained in our Section entitled “Deformation profiles close to a mismatched protein”), or equivalently . As far as the second boundary condition at is concerned, we will treat explicitly two different cases, which correspond respectively to a fixed slope and to a free slope in , as explained in the main text of the article.

Fixed slope.

In the case where the boundary conditions in are(48)which corresponds to a strong hydrophobic coupling and a fixed slope at , we obtain:(49)where(50)Note that and are either both real or complex conjugate (like ), which ensures that the solution Eq. 46 is real.

Free slope.

An alternative choice of boundary conditions in is(51)to first order again. The first of these conditions corresponds to a strong hydrophobic coupling, as before. The second one arises from minimizing the total free energy of the system without further constraints. It corresponds to the case where the slope at is free to adjust itself to yield the smallest deformation energy. With these “free-slope” boundary conditions, we obtain:(52)which are, again, either both real or complex conjugate.

Let us now assume that , as in the main text of this article. In order to understand the impact of (i.e., of ) on in the free-slope case, let us express as a function of , , and of the bulk constants , and , whose values can be extracted from the fluctuation spectra in simulations. Using Eq. 47, the relation , which can be derived from Eqs. 7 and 5, and the relation , which stems from Eqs. 4 and 7, we obtain:(53)For fixed values of , , , and , the constants can be viewed simply as functions of and : let us denote them by . The following relation holds for all and :(54)with(55)Hence, in the framework of a model that assumes , the effect of a nonvanishing on the equilibrium membrane thickness profile would be that is replaced by a renormalized spontaneous curvature , which depends linearly on . At vanishing applied tension (in which case, ), and neglecting the difference between and , we obtain Eq. 16.

3.2 Deformation energy.

Let us now calculate the deformation energy of the membrane due to the presence of the mismatched protein. For the equilibrium shape of the membrane, which is solution to the Euler-Lagrange equation Eq. 45, we are left only with boundary terms at the inclusion edge in (no other boundary terms contribute, since the deformation caused by the presence of the mismatched channel vanishes sufficiently far away from it). We can write(56)where . We have used the expression of the Gaussian curvature for small deformations in a system with cylindrical symmetry: . To express the deformation energy explicitly, one has to use the boundary conditions in .

Fixed slope.

For the boundary conditions in Eq. 48, corresponding to a fixed slope in , using Eqs. 46, 47 and 49, we can rewrite the deformation energy of the membrane in Eq. 56 as(57)This expression shows that is a second-order polynomial in and .

Spring constant for .

In the particular case where the fixed slope vanishes, Eq. 57 becomes(58)where the effective spring constant reads(59)

Dependence on applied tension.

Since , Eq. 57 shows that is a second-order polynomial in the applied tension . (In our model, features a contribution coming from , see Eq. 4. However, as mentioned in the main text, the dependence of on is negligible in practice, and we thus disregard it: in this framework, and do not depend on .) We can write(60)with(61)(62)where is the effective spring constant expressed in Eq. 59. Note that and do not appear in the coefficients and , and that and are only present in .

Free slope.

For the boundary conditions in Eq. 51, corresponding to a free slope in , using Eqs. 46, 47 and 52, we can rewrite the deformation energy of the membrane (see Eq. 56) as(63)This expression shows that is a second-order polynomial in .

Spring constant.

Eq. 63 can be expressed as(64)where the effective spring constant reads(65)while denotes the value of that minimizes , and is the minimum of , obtained for . Note that both and are nonzero if (see Eq. 63), due to the spontaneous curvature of each monolayer. The effect of monolayer spontaneous curvature was disregarded in Ref. [18], [22], which explains why Eq. 64 differs from the standard expression [22].

Dependence on applied tension.

Since , Eq. 63 shows that is a second-order polynomial in the applied tension (neglecting the -dependence of as explained in the main text). Thus, we can write(66)with(67)(68)where is the effective spring constant expressed in Eq. 65.

4 Estimating

Let us start from the free energy per molecule in monolayerexpressed in Eq. 1. All the quantities involved in this expression are defined on the hydrophilic-hydrophobic interface of the monolayer.

Let us consider a surface parallel to , and let us call the algebraic distance from to . To second order in the small dimensionless variables and , where and denote the local principal curvatures of the monolayer (recall that and ), geometry gives [19]:(69)(70)(71)Hence, we can rewrite using variables defined on , to second order:(72)where we have neglected terms containing derivatives of order higher than two.

If is the neutral surface of the monolayer [19], by definition, the curvature and the area variations are decoupled, which entails , where denotes the algebraic distance from the neutral surface to the hydrophilic-hydrophobic interface of the monolayer. Thus, given that , , and (see Methods, Sec. 1.1), we obtain(73)

Acknowledgments

We thank Grace Brannigan for sharing with us some data corresponding to the simulations described in Ref. [15]. We thank Mark Goulian for sharing with us a notebook containing some of the original calculations of Ref. [18], and for email correspondence. We also acknowledge critical reading of our manuscript by Florent Bories.

Author Contributions

Performed the experiments: AFB DC JBF. Analyzed the data: AFB DC JBF. Contributed reagents/materials/analysis tools: AFB DC JBF. Wrote the paper: AFB DC JBF.

References

  1. 1. Mouritsen OG (2005) Life { as a matter of fat: the emerging science of lipidomics. Springer, Berlin.
  2. 2. Sackmann E (1984) Physical basis for trigger processes and membrane structures. In: Chapman D, editor, Biological Membranes, London: Academic Press, volume 5. pp. 105–143.
  3. 3. Helfrich W (1973) Elastic properties of lipid bilayers - Theory and possible experiments. Zeitschrift für Naturforschung C – Journal of Biosciences 28: 693–703.
  4. 4. Owicki JC, Springgate MW, McConnell H (1978) Theoretical study of protein-lipid interactions in bilayer membranes. Proc Natl Acad Sci USA 75: 1616–1619.
  5. 5. Owicki JC, McConnell H (1979) Theory of protein-lipid and protein-protein interactions in bilayer membranes. Proc Natl Acad Sci USA 76: 4750–4754.
  6. 6. Mouritsen OG, Bloom M (1984) Mattress model of lipid-protein interactions in membranes. Biophys J 46: 141–153.
  7. 7. Huang HW (1986) Deformation free energy of bilayer membrane and its effect on gramicidin channel lifetime. Biophys J 50: 1061–1070.
  8. 8. Kelkar DA, Chattopadhyay A (2007) The gramicidin ion channel: A model membrane protein. Biochim Biophys Acta – Biomembr 1768: 2011–2025.
  9. 9. O'Connell AM, Koeppe RE II, Andersen OS (1990) Kinetics of gramicidin channel formation in lipid bilayers: transmembrane monomer association. Science 250: 1256–1259.
  10. 10. Lundbaek JA, Koeppe RE II, Andersen OS (2010) Amphiphile regulation of ion channel function by changes in the bilayer spring constant. Proc Natl Acad Sci USA 107: 15427–15430.
  11. 11. Helfrich P, Jakobsson E (1990) Calculation of deformation energies and conformations in lipid membranes containing gramicidin channels. Biophys J 57: 1075–1084.
  12. 12. Dan N, Pincus P, Safran SA (1993) Membrane-induced interactions between inclusions. Langmuir 9: 2768–2771.
  13. 13. Aranda-Espinoza H, Berman A, Dan N, Pincus P, Safran S (1996) Interaction between inclusions embedded in membranes. Biophys J 71: 648.
  14. 14. Brannigan G, Brown FLH (2006) A consistent model for thermal uctuations and protein-induced deformations in lipid bilayers. Biophys J 90: 1501–1520.
  15. 15. Brannigan G, Brown FLH (2007) Contributions of gaussian curvature and nonconstant lipid volume to protein deformation of lipid bilayers. Biophys J 92: 864–876.
  16. 16. West B, Brown FLH, Schmid F (2009) Membrane-protein interactions in a generic coarse-grained model for lipid bilayers. Biophys J 96: 101–115.
  17. 17. Elliott JR, Needham D, Dilger JP, Haydon DA (1983) The effects of bilayer thickness and tension on gramicidin single-channel lifetime. Biochim Biophys Acta – Biomembr 735: 95–103.
  18. 18. Goulian M, Mesquita ON, Fygenson DK, Nielsen C, Andersen OS, et al. (1998) Gramicidin channel kinetics under tension. Biophys J 74: 328–337.
  19. 19. Safran SA (1994) Statistical thermodynamics of surfaces, interfaces and membranes. Addison-Wesley
  20. 20. Bitbol AF, Peliti L, Fournier JB (2011) Membrane stress tensor in the presence of lipid density and composition inhomogeneities. Eur Phys J E 34: 53.
  21. 21. Fournier JB (1999) Microscopic membrane elasticity and interactions among membrane inclusions: interplay between the shape, dilation, tilt and tilt-difference modes. Eur Phys J B 11: 261–272.
  22. 22. Nielsen C, Goulian M, Andersen OS (1998) Energetics of inclusion-induced bilayer deformations. Biophys J 74: 1966–1983.
  23. 23. Watson MC, Penev ES, Welch PM, Brown FLH (2011) Thermal uctuations in shape, thickness, and molecular orientation in lipid bilayers. J Chem Phys 135: 244701.
  24. 24. Israelachvili JN (1992) Intermolecular and surface forces, Second edition. Academic Press.
  25. 25. Sharp KA, Nicholls A, Fine RF, Honig B (1991) Reconciling the magnitude of the microscopic and macroscopic hydrophobic effects. Science 252: 106–109.
  26. 26. May S, Ben-Shaul A (1999) Molecular theory of lipid-protein interaction and the Lα-H transition. Biophys J 76: 751–767.
  27. 27. Hladky SB, Gruen DWR (1982) Thickness uctuations in black lipid membranes. Biophys J 38: 251–258.
  28. 28. Rawicz W, Olbrich KC, McIntosh T, Needham D, Evans E (2000) Effect of chain length and unsaturation on elasticity of lipid bilayers. Biophys J 79: 328–339.
  29. 29. May ER, Narang A, Kopelevich DI (2007) Role of molecular tilt in thermal uctuations of lipid membranes. Phys Rev E 76: 021913.
  30. 30. May ER, Narang A, Kopelevich DI (2007) Molecular modeling of key elastic properties for inhomogeneous lipid bilayers. Mol Simulat 33: 787–797.
  31. 31. Watson MC, Brandt EG, Welch PM, Brown FLH (2012) Determining biomembrane bending rigidities from simulations of modest size. Phys Rev Lett 109: 028102.
  32. 32. Lindahl E, Edholm O (2000) Mesoscopic undulations and thickness uctuations in lipid bilayers from molecular dynamics simulations. Biophys J 79: 426–433.
  33. 33. Marrink SJ, Mark AE (2001) Effect of undulations on surface tension in simulated bilayers. J Phys Chem B 105: 6122–6127.
  34. 34. Lundbaek JA, Andersen OS (1999) Spring constants for channel-induced lipid bilayer deformations estimates using gramicidin channels. Biophys J 76: 889–895.
  35. 35. Kim T, Lee KI, Morris P, Pastor RW, Andersen OS, et al. (2012) Inuence of hydrophobic mismatch on structures and dynamics of gramicidin A and lipid bilayers. Biophys J 102: 1551–1560.
  36. 36. Kolb HA, Bamberg E (1977) Inuence of membrane thickness and ion concentration on the properties of the gramicidin A channel: autocorrelation, spectral power density, relaxation and singlechannel studies. Biochim Biophys Acta – Biomembr 464: 127–141.
  37. 37. Hwang TC, Koeppe RE II, Andersen OS (2003) Genistein can modulate channel function by a phosphorylation-independent mechanism: importance of hydrophobic mismatch and bilayer mechanics. Biochemistry 42: 13646–13658.
  38. 38. Chung H, Caffrey M (1994) The curvature elastic-energy function of the lipid-water cubic mesophase. Nature 368: 224–226.
  39. 39. Chung H, Caffrey M (1994) The neutral area surface of the cubic mesophase: location and properties. Biophys J 66: 377–381.
  40. 40. Templer RH (1995) On the area neutral surface of inverse bicontinuous cubic phases of lyotropic liquid crystals. Langmuir 11: 334–340.
  41. 41. Templer RH, Turner DC, Harper P, Seddon JM (1995) Corrections to some models of the curvature elastic energy of inverse bicontinuous cubic phases. J Phys II France 5: 1053–1065.
  42. 42. Vacklin H, Khoo BJ, Madan KH, Seddon JM, Templer RH (2000) The bending elasticity of 1-monoolein upon relief of packing stress. Langmuir 16: 4741–4748.
  43. 43. Szule JA, Fuller NL, Rand RP (2002) The effects of acyl chain length and saturation of diacylglycerols and phosphatidylcholines on membrane monolayer curvature. Biophys J 83: 977–984.
  44. 44. Leikin S, Kozlov MM, Fuller NL, Rand RP (1996) Measured effects of diacylglycerol on structural and elastic properties of phospholipid membranes. Biophys J 71: 2623–2632.