Bilayer Elasticity at the Nanoscale: The Need for New Terms

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][5][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][12][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][15][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][15][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 z and {, in contact with a reservoir of lipids with chemical potential m. We write the effective Hamilto-nian per molecule in monolayer + as where S + is the area per lipid, while H + is the local mean curvature of the monolayer, and K + is its local Gaussian curvature (denoting by c + 1 and c + 2 the local principal curvatures [19] of the monolayer, we have H +~( c + 1 zc + 2 )=2 and K +~c+ 1 c + 2 ). All these quantities are defined on the hydrophilichydrophobic interface of each monolayer. Eq. 1 corresponds to an expansion of f + m for small deformations around the equilibrium state where the membrane is flat and where each lipid has its equilibrium area S 0 . Any constant term in the free energy per lipid is included in a redefinition of the chemical potential m. From now on, we will consider small deformations of an infinite flat membrane and we will work in the Monge gauge, so H +^+2 h + =2 and K +^L2 x h + L 2 y h + {(L x L y h + ) 2~d et (L i L j h + ), where z~h + (x,y) represents the height of the hydrophilichydrophobic interface of each monolayer with respect to a reference plane (x,y). The upper monolayer is labeled by z and the lower one by {. Many constants involved in Eq. 1 can be related to the constitutive constants of a monolayer: f '' 0 S 0~Ka =2 is the compressibility modulus of the monolayer, f 2 =(2S 0 )~k 0 =2 is its bending rigidity, f K =S 0~ k k=2 is its Gaussian bending rigidity, f 1 =f 2~c0 is its spontaneous (total) curvature, and f ' 1 =f 2~c ' 0 is the modification of the spontaneous (total) curvature due to area variations (see Methods, Sec. 1.1).
In the case where a~b~f~0, Eq. 1 is equivalent to the model of Ref. [19], which is the basis of that developed in Refs. [12][13][14][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 a, b and f cannot be neglected with respect to others. We focus on the influence of a, for which we provide a physical interpretation, and we will set b~f~0 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 a, b and f are all included.
The effective Hamiltonian of a bilayer membrane patch with projected area A p at chemical potential m 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 h~(h z zh { )=2 of the bilayer, the sum u of the excess hydrophobic thicknesses of the two monolayers, each being measured along the normal to the monolayer hydrophilichydrophobic interface (see Fig. 1 and Eqs. [26][27][28][29], and the difference d 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.) In the present work, we are not interested in the degree of freedom d, 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 d out, which amounts to minimizing f with respect to d since our theory is Gaussian. The resulting effective Hamiltonian, which involves h and u, is given by Eq. 32 in Sec. 1.2 of our Methods part. In this effective Hamiltonian, the variables h and u are decoupled, and the part depending on h 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 h (see Methods, Sec. 1.3).
Here, we focus on variations of the membrane thickness, i.e., on the variable u. We thus restrict ourselves to the case where the average shape h of the membrane is flat (see Fig. 1B). In this case, we obtain, from Eq. 32: 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 u z zd 0 =2 (resp. u { zd 0 =2). The height of monolayer + along z is denoted by h + . The average membrane shape, h~(h z zh { )=2, 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 d 0 , 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 (u z~u{~u =2). Hence, the average shape h is constant, and confounded with the midlayer of the membrane. Although u + is defined along the normal to the monolayer hydrophilic-hydrophobic interface, the boundary condition at the inclusion edge, i.e., in r~r 0 , simply reads u(r 0 )~u 0~' {d 0 to first order (see main text, Section entitled ''Deformation profiles close to a mismatched protein''). doi:10.1371/journal.pone.0048306.g001 In the case where b~f~0, on which the body of this article focuses, the various constants introduced in Eq. 2 read: K'' a~k In these equations, d 0 denotes the equilibrium hydrophobic thickness of the bilayer membrane, s plays the part of an externally applied tension (see Methods, Sec. 2), K a is the compressibility modulus of the membrane, k k is its Gaussian bending rigidity, k 0 is the bending rigidity of a symmetric membrane such that d~0, c 0 is the spontaneous (total) curvature of a monolayer, and c' 0 is the modification of this spontaneous curvature due to area variations. In addition, we have introduced k' a~4 a S 0 =d 2 0 , which has the dimension of a surface tension, like K a . Note that the last three terms in Eq. 2 are boundary terms.
In Sec. 1.2 of our Methods part, the expressions of K' a , K'' a , A 1 and A 2 are provided in the more general case where b and f 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 K a w0, K'' a w0, and

Discussion
Comparison with existing models Our model Eq. 2 has a form similar to that of the models developed in Refs. [12][13][14][15]. However, it differs from these previous models on several points. First, our definition of u is slightly different. Second, we have included the effect of an applied tension s. 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 u. In the present work, the variable u, which is the relevant one to study membrane thickness deforma-tions, 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. [26][27][28][29] in the Methods section). This definition of u has the advantage of being independent of deformations of the average membrane shape h.
The excess thickness variable used in Refs. [7,[12][13][14][15]18,22,23] reads in our notations: Using Eqs. 9 and 25, and working to second order, we obtain which shows that there is a second-order difference between 2 u u and our variable u. Consequently, the difference between the definition used in the previous works and ours regards only the term linear in u, 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 u or 2 u u. Our definition of u is the right one for rigorously taking tension into account, because it is independent of deformations of the average membrane shape h: the energy stored in the variable u only comes from thickness variations. (The variable u u of Refs. [7,[12][13][14][15]18,22,23] corresponds to the difference between the bilayer hydrophobic thickness projected along z and the non-projected equilibrium hydrophobic bilayer thickness (see Eq. 9), so it is not independent of h. The second-order difference between 2 u u and u, which is shown in Eq. 10, arises from this difference in projection between actual thicknesses and equilibrium thicknesses within the definition of u u.) On tension. First of all, existing models [7,[12][13][14][15]18,22] were constructed at zero applied tension, which means s~0 in Eq. 2. To our knowledge, our work is the first where the coefficient of the term linear in u 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 K' a . In our model, the constant K' a 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 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, k' a , arises from a, i.e., from the term in (+S) 2 introduced in Eq. 1. This term was not included in Refs. [12][13][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 c 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 (u z~u{~u =2) with flat average shape, the surface of the hydrocarbon-water interface is increased by a factor ½1z(+u) 2 =8 for each monolayer (see Fig. 2). Thus, to second order, the associated energetic cost per unit projected area is c(+u) 2 =4. Note that other physical effects, e.g., the elasticity of the chains, may yield contributions to the term in (+S) 2 . However, if we restrict to the simple term arising from interfacial tension, we obtain Finally, the third contribution, s=4, 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 k' a in Eq. 12, we expect s=4%k' a .
In the seminal article Ref. [7], where the membrane model was constructed by analogy with liquid crystals, a term in (+u) 2 , 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 (+u) 2 .
However, our construction of the membrane effective Hamiltonian shows that the microscopic tension involved through k' a 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 shortwavelength 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 (+u) 2 is negligible if where we have used our own notations of the prefactors of the terms in (+u) 2 , u 2 and (+ 2 u) 2 . In the case of DOPC, taking K'' a~k =4 and using the values of the membrane constants [28], this condition becomes K' a %28 mN=m. While this is well verified if K' a 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 K' a , in particular the one arising from interfacial tension. Besides, in Sec. 1.2 of our Methods part, we show that b is also involved in K' a , 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 k' a , and comparing the relative weight of the different contributions to K' a .
On the value of K'' a . We have obtained K'' a~k0 =4 (see Eq. 5), where k 0 is the bending rigidity of a symmetric membrane such that d~0. The elastic constant k 0 is related to the bending rigidity k of the Helfrich model (see Methods, Sec. 1.3) through The difference between k 0 and k arises from integrating out d (see Methods, Sec. 1.2). In the previous models, this procedure was not carried out, as one focused directly on the symmetric case d~0.
All previous models thus made the approximation K'' a~k =4 [7,[12][13][14]18,22]. In addition, in Sec. 1.2 of our Methods part, we show that f is also involved in K'' a , 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 Bilayer Elasticity at the Nanoscale PLOS ONE | www.plosone.org first two boundary terms are the same as in Refs. [12][13][14]. However, even at vanishing applied tension, we have K' a ={2 A 2 , contrary to the previous models [14], due to the presence of k' a . We have also accounted for the Gaussian bending rigidity k k, as in Ref. [15]: it yields the third boundary term.
Again, the situation is more complex when b is included, as the expressions of A 1 and A 2 then feature extra terms linear in b (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][30][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 s~0. 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 u).
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 k, while those regarding the thickness of the membrane give access to K a , K' a and K'' a .
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 K'' a~k =4, in order to include the possible effect of the difference between k and k 0 (see Eq. 33), and of f (see Eq. 35). Our results were similar to those obtained in Refs. [14,16] assuming that K'' a~k =4, and we obtained no systematic significant difference between K'' a and k=4, which means that the corrections to K'' a predicted by our model are negligible in these simulations. This gives a justification for focusing only on the correction to K' a , as we do in this article. Besides, we obtained K' a v0 from all the fits, as reported in Refs. [14,16], and we checked that all the values obtained for K' a comply with the stability condition Eq. 8.
Deformation profiles close to a mismatched protein. In Refs. [14][15][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 k' a .
Let us denote the radius of the protein by r 0 , and its hydrophobic length by ': the mismatch originates from the difference between ' and the equilibrium hydrophobic thickness d 0 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 r~r 0 . There is a consensus on the assumption of strong hydrophobic coupling u(r 0 )~u 0~' {d 0 , 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 u, the condition u(r 0 )~u 0~' {d 0 is valid to first order, while it is exactly valid with the definition of Refs. [7,[12][13][14][15]18,22,23] (see Eqs. 9, 10). This difference arises from the fact that our u is not projected along z (see Fig. 1), which makes it fully independent of h. 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 r 0 (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 r~r 0 [12][13][14][15][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 f s favoring a slope s 0 in r 0 : where k s is an effective rigidity, while u' denotes the derivative of the membrane thickness profile u with respect to the radial coordinate r. Two a priori unknown parameters, k s and s 0 , are associated with this effective potential. The ''free-slope'' boundary condition (also called ''natural'' boundary condition [12,14]) is recovered in the limit k s ?0, which is appropriate if f s is negligible with respect to the energetic contributions in f . Conversely, if k s ??, the protein locally imposes the fixed slope s 0 . 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][15][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][15][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][15][16] is very similar to ours, except that k' a~0 . 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 bỹ c c 0 , was found to depend linearly on the hydrophobic mismatch u 0 [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 k' a . We show in Sec. 3.1 of our Methods part that the quantityc then plays the part of the renormalized spontaneous curvature of Ref. [16] in the equilibrium profile. This quantity is linear in u 0 : our model, and more precisely the presence of a nonvanishing k' a , thus provides an appealing explanation for the linear dependence observed in Ref. [16]. Using a linear fit of the data of Ref. [16] (see Fig. 4), together with Eq. 16 and the value k~2:3|10 {20 J extracted from the spectra in Ref. [16], we obtain k' a~1 3 mN=m.
It is interesting to compare this value to K' a , which is obtained from the fluctuation spectra in Ref. [16]: K' a~{ 9:2mN=m. This shows that the contribution of k' a to K' a is important. Besides, we may now estimate the contribution to K' a that stems from the monolayer spontaneous curvature (see Eq. 4): Using values from the fluctuation spectra in Ref. [16], this yields j&{6A 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 j and c' 0 ).
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 r 0 . 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 A + (0,c c 0 ) (see Eq. 54), usingc c 0 as our only fitting parameter, in the spirit of Ref. [16]. We found that c c 0 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 u 0 , we have averagedc c 0 over the seven results corresponding to the different inclusion radii. The line joining these two average values ofc c 0 as a function of u 0 is plotted in Fig. 5B. Using Eq. 16 and the value k~1:4|10 {19 J [14,15], the slope of this line yields k' a~3 6 mN=m: the order of magnitude of this value is the same as the one obtained from the data of Ref. [16].
Again, we can compare this value to K' a , which is obtained from the fluctuation spectra in Refs. [14,15]: K' a~{ 11:9mN=m. Hence, the contribution of k' a to K' a is important here too. We also obtain {k 0 (c 0 {c' 0 S 0 )=d 0~K ' a {k' a~{ 48 mN=m, and j&{3A.
In Ref. [15], the shortcomings of the model that disregards k' a are explained by the local variation of the volume per lipid close to  [16] in the vicinity of a mismatched inclusion with hydrophobic thickness '^2:4nm and radius r 0~9 A, with center in r~0, as a function of the radial coordinate r. The equilibrium membrane hydrophobic thickness is d 0^3 :6nm. 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. 46-53 with k' a~0 , taking u 0 and the (renormalized) spontaneous curvaturec c 0 as fitting parameters, the other constants being known from the fluctuation spectra. doi:10.1371/journal.pone.0048306.g003 Figure 4. Renormalized spontaneous curvaturec c 0 as a function of the hydrophobic mismatch u 0 . Data from Ref. [16], which presents fits of simulation results for inclusions with three different hydrophobic thicknesses. Line: linear fit, with slope (0:56+0:02)nm {2 . Note that ourc c 0 corresponds to twice that in Table 2 where v 0 is the bulk equilibrium volume per lipid, while v 0 zv(r 0 ) denotes the volume per lipid in r~r 0 . However, the predicted linear dependence ofc c 0 in v(r 0 )=v 0 is not obvious: in Fig. 6, we rather see two groups of points (one for each value of u 0 ) than a linear law. In other words, the data of Ref. [15] is more consistent with a value ofc c 0 that depends only on u 0 and not on v (or r 0 ), 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 ofc c 0 versus u 0 [16].
To sum up, our model can explain the dependence ofc c 0 in u 0 observed in the numerical results of Refs. [15,16] as a consequence of the presence of k' a . 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 k' a 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 Nterminus 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 updown 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][12][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 r~r 0 . As discussed in the previous section, we will assume strong hydrophobic coupling, i.e., u(r 0 )~u 0~' {d 0 , 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 F associated with the deformation can be expressed as F~Hu 2 0 , where H is the effective spring constant, while u 0 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].   [15], and the values of v(r 0 )=v 0 are directly taken from Ref. [15]. doi:10.1371/journal.pone.0048306.g006 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 H 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 H exp~1 15+10 mN:m {1 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 u 0 . The prefactor of u 2 0 in the deformation energy F 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, k' a is involved in H, through K' a . Our aim will be to find out which value of k' a gives the best agreement with the experimental value of H.
Using Eqs. 4, 5 and 7, and neglecting the difference between k and k 0 , Eqs. 59 and 65 show that H depends on the elastic constants k, k k and c 0 involved in the Helfrich model, on K a , on c' 0 S 0 , which corresponds to the spontaneous curvature variation with the area per lipid, on d 0 , on the radius r 0 of the gramicidin channel, and on k' a . There is, to our knowledge, no direct experimental measurement of c' 0 S 0 available, but, as shown in Sec. 4 our Methods part, we have c' 0 S 0~Ka j=k, where j 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 k and k 0 ). Hence, in order to calculate the spring constant, we need values for k, k k, c 0 , K a and j, in the precise case of monoolein membranes.
In Ref. [38], the elastic constants k, k k and c 0 were measured in a monoolein cubic mesophase, both at 25 0 C and at 35 0 C. 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 c 0 , and on the ratio k k=k [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 H II phase, where the strains are better relaxed. However, these measurements were carried out at 37 0 C, while the experiments of Ref. [17] that we wish to analyze were performed at 23 0 C. 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 K a for monoolein is provided by Ref. [27].
In Table 1, we present the results obtained for the spring constant H 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 k k [38,41]. However, k k is involved in H only in the free-slope case (see Eqs. 59 and 65): the 3% difference between the values of H 0 obtained with parameter sets 1 and 2 stems only from the difference on c 0 , while the 12% difference between H f obtained with data sets 1 and 2 contains an important contribution from k k. 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 H 0 and a 20% difference on H f . 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 (j~{0:3A [40]), the contribution of c' 0 S 0 to H 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 H f obtained assuming that k' a~0 are about three times smaller than the experimental value H exp~1 15+10 mN:m {1 (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 H.) However, H f reaches the experimental value for k' a &25 mN=m for all three parameter sets (see line 2 of Table 1). Hence, for free-slope boundary conditions, the presence of k' a , with an order of magnitude consistent with Eq. 12, improves the agreement between theory and experiment.
We may compare these values of k' a to the contribution to K' a that originates from the monolayer spontaneous curvature (see Eq. 4): {k 0 (c 0 {c' 0 S 0 )=d 0 . We estimate the value of this contribution to be between 0:26 and 1:2mN=m, 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 j seemed to be significant in the numerical simulations. In addition, the contribution of membrane tension to K' a , namely, s=4, cannot The results are given both for the free-slope boundary condition (using Eq. 65) and for the zero-slope boundary condition s~0 (using Eq. 59). All values of H and k' a are given in mN=m. Negative values of k' a are not detailed since they would yield an instability for the monolayer Hamiltonian Eq. 22 in the present framework where k'' a~0 . The different columns correspond to three different data sets for the parameters of the membrane. Set 1 corresponds to the data from [38] at 25 0 C: k~3:6|10 {20 J, c 0~{ 0:135 nm {1 , k k~8:8|10 {22 J. Set 2 takes into account the corrections on c 0 and k k in [41]: c 0~{ 0:196 nm {1 , k k~{3:6|10 {21 J. Set 3 corresponds to the data from [42]: k~1:2|10 {20 J, c 0~{ 0:503 nm {1 , and k k~{1:2|10 {21 J deduced from k k=k~{0:1 [41]. In all cases, we have taken r 0~1 nm [34], d 0~2 :46nm [39], j~{0:3A [40], K a~1 40 mN=m [27,34]. doi:10.1371/journal.pone.0048306.t001 exceed about 1 mN/m. In the case of the free-slope boundary condition, our results imply that k' a should be the dominant contribution to K' a 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 H 0 assuming that k' a~0 are in quite good agreement with the experimental value H exp~1 15+10 mN:m {1 obtained in Ref. [34] from the data of Ref. [17], for all the data sets we used (see line 3 of Table 1): hence, k' a 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 f in lipid vesicles as a function of the tension s 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: ln f~g(s)~C 0 zC 1 szC 2 s 2 : ð18Þ Given that ln f 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 s, since the membrane hydrophobic thickness depends on s (see Eq. 43 in Sec. 2 of our Methods part). Expressing the deformation energy F of the membrane due to the presence of the dimer gives a theoretical expression for the s-dependent part of ln f . In our model, k' a features a contribution coming from s (see Eq. 4). However, this term is negligible, given that s=4% ffiffiffiffiffiffiffiffiffiffiffiffi K a K'' a p =d 0 (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 F on s, in agreement with the form of Eq. 18. Comparing the experimental values of C 1 and C 2 to those predicted by theory provides a test for theoretical models [18]. (Note that, if the sdependent contribution to K' a is included, the expression of the sdependent part of ln f is no longer simply quadratic in s. However, we explicitly verified that including this contribution yields a negligible change to the relation between s and ln f , for realistic values of the parameters).
Since the coefficients C 1 and C 2 arise from membrane elasticity, they are common to all the vesicles studied in Ref. [18], which have the same lipid composition. Conversely, the baseline C 0 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 where the index j runs over all the experimental points, with fitting parameters C 1 ,C 2 ,C k 0 ,k~1, . . . ,12. The baseline C k 0 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 In Ref. [18], the authors used published values of the material constants to calculate C 1 and C 2 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 s (s~0:3). However, we need to raise the following points: 1. There was a mistake in their implementation of the formula of Ref. [22] giving C 1 and C 2 as a function of the material constants. More precisely, we found that a factor of 2 was missing in the expression of C 1 and a factor of 4 was missing in that of C 2 in the implementation of the formula of Ref. [22]. This was confirmed by Mark Goulian (private communication).
The actual values of C 1 and C 2 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: d 0~2 :7nm [18], K a~2 65mN=m, k~8:5|10 {20 J [28], c 0~{ 0:132nm {1 [43], and the dimensions of a gramicidin channel: r 0~1 nm, ''~'zd~2:3nm [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).
A somewhat better agreement with the experimental data is obtained when taking s~0 instead of s~0:3 for the fixed slope (see Figs. 7 and 8, and Table 2). However, the downward inflection of the experimental curves at high s is not adequately described for any value of s. In fact, C 2 is independent of s, 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 F in the framework of our model, both for the fixed-slope boundary condition and for the free-slope boundary condition. The resulting expressions of C 1 and C 2 are given by Eqs. 61, 62, 67 and 68. In order to see which values of k' a and which boundary conditions give the best agreement with the experiments of Ref. [18], we present a plot of the goodness-of-fit function x 2 (see Eq. 19) in a (C 1 ,C 2 ) graph in Fig. 8. On this graph, we have plotted the trajectories obtained from our model in the (C 1 ,C 2 ) plane when varying k' a , for s~0, for s~0:3 (as in Ref. [18]), and for the free-slope boundary condition.
In order to obtain numerical values of C 1 and C 2 from Eqs. 61, 62, 67 and 68, we used the above-mentioned parameter values, and the estimate k k~{0:8k [16]. Finally, we estimated c' 0 S 0 through the relation c' 0 S 0~Ka j=k (see Eq. 73 in Sec. 4 of our Methods part). For this, the algebraic distance j 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 j&{0:5A. Here again, the neutral surface is close to the hydrophilic-hydrophobic interface. For the sake of simplicity, we took c' 0~0 , and we checked that the results were not significantly different when taking j~{0:5A.
The ingredient in our model that can change significantly the results is k' a (Note that the values of C 1 and C 2 corresponding to k' a~0 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 C 1 can be explained by our model. In addition, the values of k' a that minimize x 2 , i.e., that give the best agreement with the experimental data of Ref. [18], are between 0 and 50 mN=m, depending on the boundary condition chosen, as shown in Table 3. This range of values of k' a is reasonable.
For the free-slope boundary condition, the best agreement with the experimental results is obtained for k' a &40 mN=m (see Table 3 and Fig. 8). The order of magnitude is the one expected from k' a &c=2.
Let us now discuss the results obtained for the fixed-slope boundary condition, which is used in Ref. [18]. For a fixed slope s~0, the best agreement with the results of Ref. [18] analyzed with the complete quadratic fit is obtained for k' a~0 . Conversely, for s~0:3, the best agreement is obtained for k' a &40 mN=m, 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 s that is chosen.
In all cases, the absolute values of C 2 we obtain remain much smaller than the one that matches best the experimental results, which is C 2~{ 90:0|10 {3 (mN=m) {2 (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 C 2 . Our model predicts that C 2 is proportional to the effective spring constant H 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 H but not with those of C 2 . This disagreement on C 2 could come either from a shortcoming of the model or from an undetected systematic error in the experimental data. The importance of C 2 is largest at highest tensions, as it is C 2 which gives the curve its concavity, and it should be noted that the maximum applied tension s is around 4:5mN=m in Ref. [18], which is comparable to the rupture threshold of 3{10 mN=m [28]. The membrane properties may be affected at such high tensions in a way that is no longer well described by standard  Table 2. doi:10.1371/journal.pone.0048306.g007 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 C 2 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 C 2~0 ), keeping only the points corresponding to sv2mN=m: this yields C 1~( 0:62+0:05) (mN=m) {1 (see Fig. 9). In Table 4, we list, for different boundary conditions, the value of k' a which gives C 1~0 :62 (mN=m) {1 , and the value of C 2 obtained from our model for this k' a . 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 k' a 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 100mN=m. Again, these values depend a lot on s for fixedslope boundary-conditions. (For instance, the slope s~{0:17 is consistent with k' a~0 (see Table 4). However, there is no a priori reason for assuming that k' a~0 .) Again, we may compare our estimates of k' a (see Tables 3 and 4) to the term {k 0 (c 0 {c' 0 S 0 )=d 0 , which also contributes to k' a : here, {k 0 (c 0 {c' 0 S 0 )=d 0~{ 0:76 mN=m. 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 Figure 8. Comparison between the experimental values of C 1 and C 2 and those obtained from different models. Colorscale: goodnessof-fit function x 2 (see Eq. 19) for the data of Ref. [18], as a function of the fitting parameters C 1 and C 2 . White diamond: values of C 1 and C 2 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 (C 1 ,C 2 ) plane when varying k' a . Red: free slope; green: s~0, black: s~0:3. These three curves start by a white dot at k' a~0 , and k' a increases rightwards along these curves. The rightmost white dot (k' a~0 , s~0) 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). doi:10.1371/journal.pone.0048306.g008 Table 2. Values of C 1 and C 2 obtained from the model of Ref. [22] and from our model with k' a~0 .

Model
Ref. [22] Ref. [22] Ref. [ The results are presented both for the fixed-slope boundary condition (see Eqs. 61 and 62), with slopes 0 and 0:3, and for the free-slope boundary condition (see Eqs. 67 and 68). The corresponding values of x 2 are also given. Recall that the best quadratic fit to the data of Ref. [18]  monolayer and its hydrophilic-hydrophobic interface are very close, while j 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 k' a should be the dominant contribution to K' a 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 k' a , which is of the order of magnitude expected if k' a~c =2. Conversely, for the fixed-slope boundary condition, the value of k' a is coupled to that of the slope s.

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 k' a in the range 13{60 mN=m. Such values are consistent with the idea that this term involves a significant contribution of the interfacial tension c between water and the hydrocarbon-like hydrophobic part of the membrane. Indeed, this contribution should yield k' a~c =2&25 mN=m.
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 Table 3. Values of k' a , C 1 and C 2 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).  Figure 9. Formation rate f of gramicidin channels as a function of the applied tension s, analyzed with a linear model, for sv2mN=m. Diamonds: experimental data retrieved from Fig. 6b of Ref. [18], after subtraction of the baselines C k 0 (which are different from those of Fig. 7 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: where K n is the n th -order modified Bessel function of the second kind, and which are either both real or complex conjugate. The integration constants A + are determined by the boundary conditions at r~r 0 . 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 u(r 0 )~u 0~' {d 0 (to first order, as explained in our Section entitled ''Deformation profiles close to a mismatched protein''), or equivalentlyũ u(r 0 )~ũ u 0~' {d 0 1{s=K a ð Þ . As far as the second boundary condition at r~r 0 is concerned, we will treat explicitly two different cases, which correspond respectively to a fixed slope and to a free slope in r 0 , as explained in the main text of the article.
Fixed slope. In the case where the boundary conditions in r~r 0 areũ u(r 0 )~ũ u 0~' {d 0 1{ s K a u u'(r 0 )~s which corresponds to a strong hydrophobic coupling and a fixed slope s at r~r 0 , we obtain: