Gramicidin A Channel Formation Induces Local Lipid Redistribution II: A 3D Continuum Elastic Model

Article

Gramicidin A Channel Formation Induces Local Lipid Redistribution II: A 3D Continuum Elastic Model Alexander J. Sodt,1,* Andrew H. Beaven,2 Olaf S. Andersen,3 Wonpil Im,4 and Richard W. Pastor5 1

Eunice Kennedy Shriver National Institute of Child Health and Human Development, National Institutes of Health, Bethesda, Maryland; Department of Chemistry, The University of Kansas, Lawrence, Kansas; 3Department of Physiology and Biophysics, Weill Cornell Medical College, New York, New York; 4Department of Biological Sciences and Bioengineering Program, Lehigh University, Bethlehem, Pennsylvania; and 5Laboratory of Computational Biology, National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, Maryland 2

ABSTRACT To change conformation, a protein must deform the surrounding bilayer. In this work, a three-dimensional continuum elastic model for gramicidin A in a lipid bilayer is shown to describe the sensitivity to thickness, curvature stress, and the mechanical properties of the lipid bilayer. A method is demonstrated to extract the gramicidin-lipid boundary condition from allatom simulations that can be used in the three-dimensional continuum model. The boundary condition affects the deformation dramatically, potentially much more than typical variations in the material stiffness do as lipid composition is changed. Moreover, it directly controls the sensitivity to curvature stress. The curvature stress and hydrophobic surfaces of the all-atom and continuum models are found to be in excellent agreement. The continuum model is applied to estimate the enrichment of hydrophobically matched lipids near the channel in a mixture, and the results agree with single-channel experiments and extended molecular dynamics simulations from the companion article by Beaven et al. in this issue of Biophysical Journal.

INTRODUCTION The lipid environment must adapt to the specific conformation of the protein it surrounds. The same strong hydrophobic interaction that localizes proteins to the membrane also ensures that the thickness of the protein and membrane are matched (1,2). The soft lipid leaflets compress (3) and bend (4–8), and their lipids tilt (9–11) to accommodate the protein, including its cavities and convexities (12). Thus a conformation can be resisted or promoted by its compatibility with the surrounding bilayer (12–15). The membrane itself has been successfully treated as a material, albeit one with significant vertical inhomogeneity and thus stress (16,17). This is accomplished through a continuum elastic model (CEM), described by elastic moduli (18) and local stresses, and with boundary conditions enforcing the strong hydrophobic matching as a constraint. A fundamental assumption of the elastic model of the mem-

Submitted September 27, 2016, and accepted for publication January 30, 2017. *Correspondence: [email protected] Alexander J. Sodt and Andrew H. Beaven contributed equally to this work. Editor: Scott Feller. http://dx.doi.org/10.1016/j.bpj.2017.01.035

1198 Biophysical Journal 112, 1198–1213, March 28, 2017

brane is that all the complexity of the internal conformation of the protein can be summarized in this boundary condition, to which the membrane adapts across all protein classes. In the material view, two problems emerge: how lipid composition affects the local elastic parameters (19– 22) and, once a model is determined, what force a local lipid composition exerts on protein conformation (23–26) and vice versa (27,28). These questions tend to have been pursued using a twodimensional (2D) CEM to describe the membrane deformation around a protein. The model begins with the Helfrich Hamiltonian (4,5) for (lateral) curvature energetics, then adds a term modeling the vertical (or normal) compression of the material that occurs to achieve hydrophobic matching at the protein. The model is attractive because the energy function, EL (Eq. 2 below) is parameterized by experimentally available quantities. However, at the crucial interface between the three-dimensional (3D) protein and the 2D bilayer, a choice must be made for the connection. Without independent guidance, alternate (29,30), but conceptually reasonable, boundary conditions can quadruple the deformation energy for a particular mismatch.

Modeling the Gramicidin-Lipid Boundary

In this work, the 2D-CEM (31) is expanded to include the deformation of the bilayer not just at the surface as in the 2D-CEM, but throughout its thickness. The model, termed the 3D-CEM, takes into account the connection of the bilayer material and inclusion considering the entire surface of the protein. The extension comes at the expense of adding additional parameters that describe the origins of the 2D experimental parameters from local 3D properties of the material. Naturally, these choices must be validated. The two quantities evaluated from the models (including the all-atom model) are directly related to the two energetics terms at play: curvature and compression. First, the curvature stress of the leaflet is available explicitly from simulation (32) and implicitly from experiments on the response of the protein to lipid-compositional curvature stress (13,23,33). Second, the lipid redistribution around the protein is available explicitly from simulation and implicitly from experiments on the response of proteins to mixed-lipid environments. A basic tool for this and similar studies is the use of gramicidin A (gA) as a ‘‘force probe’’ (34,35). The transmembrane channel is formed from two monomers on opposite leaflets. The activation energy for dissociation of the channel is modulated by the energy of the material deformation around the complete channel, and this directly impacts the channel lifetime. This system can thus be used to validate a material model of the bilayer usable for more complex systems. In this work, the model of gA has rotational symmetry. This is reasonable (though not strictly necessary) and dramatically simplifies the numerical simulations. This is of course not the general case for proteins. In (36), for example, large membrane deformations are identified near adjacent polar/apolar residues, where hydrophobic matching cannot be completely resolved. The analysis of the bilayer deformation energy through single-channel-lifetime experiments is discussed in Article I by Beaven et al. in this issue of Biophysical Journal. For clarity, some material is repeated here, but now from the perspective of inferring the peptide-lipid boundary condition. Developing and validating bilayer continuum modeling from single-channel-lifetime experiments (without molecular simulation) relies on analysis of the difference in bilayer deformation energy between the transition state (TS) energy and the dimer energy. Essentially, the deformation energy, E, is a quadratic function of the (per-leaflet) mismatch at the boundary, u0 , given by (34)
2 2 E ¼ H 2h0  2hp ¼ Hð2u0 Þ ;

(1)

where h0 is the equilibrium thickness of the lipid leaflet, hp is half the height of the inclusion (for use in the single-leaflet model), and H is a coefficient reflecting bilayer deformation energetics. Here, all contributions to the energy, including the effect of boundary conditions, are gathered into the

phenomenological spring coefficient, H. Below, model analysis will separate terms individually. The energy difference between the TS and the dimer is then Hð4u0 d  d2 Þ, and in a transition-state model, it explains the difference in rates quantified by logðt 0 =t 1 Þ, where at u0 ¼ 0, t ¼ t 0 . To validate a model against experimental timescales requires connecting t and u0 to the energy available from the model, through the quantity H. The two adjustable parameters in the analysis are the slope of the leaflet at the boundary (which influences H) and the difference in effective hydrophobic thickness between the TS and the dimer, d. Neither of these quantities is directly available from the experiment. The choice for the boundary condition has been debated since the theory of elastic bilayer deformations was applied in 1986 (30) to analyze the observed relation between gramicidin single-channel lifetimes and lipid bilayer thickness (37). This analysis showed that the s value, the slope of the leaflet at the channel, deduced was quite different, in fact negative, from the value that would minimize the deformation energy (29); a similar conclusion was reached by Lundbaek and Andersen (38). It is in this context important that the choice of s (slope) effectively lumps all the uncertainties about the model (and the molecular details of lipid packing adjacent to the channel, which are not included in the 2D CEM) into a single adjustable parameter. As elaborated on in the discussion, inferring s from single-channel lifetimes relies on knowledge of the structure of the TS, as well as on precise values of the mechanical constants. This work is focused on the use of simulation to resolve the ambiguity of the peptide-lipid boundary. All-atom simulations are used to compute the molecular details of the boundary, as well as the curvature stress that is highly correlated with the form of the boundary. In this way, the simulation is treated as if it were an experiment, i.e., that with no adjustable parameter it has predictive power. Although the force field used here, CHARMM 36 (C36) (39), reproduces assorted mechanical properties quite well (40), the accuracy of any lipid force field is open to question (41). As is pointed out in (41), although the accuracy of the CHARMM force field for nuclear-magnetic resonance observables is good, it may result in over-fitting. How lipids at the boundary of gA adapt is difficult to validate directly by experiment, and is a source of concern. This highlights the importance of the comparison with experiment, although indirect, in Article I. In summary, although simulation information does not rely on knowledge of the TS, it does rely on the validity of the force field. Characterizing the protein-lipid interface continues to be a challenge. For example, Yoo and Cui (42) have computed spatially resolved local normal and tangential stresses around the gA channel to characterize the bilayer deformation beyond using structural information. In a review of continuum membrane modeling, Argudo et al. (43) present a method for performing hybrid continuum-atomistic simulations with a continuum membrane that provides a close

Biophysical Journal 112, 1198–1213, March 28, 2017 1199

Sodt et al.

match to atomistic modeling. This work is complementary to these studies, as it develops a technique for characterizing and validating continuum peptide-lipid boundaries against simulation. The theory in this study cross-validates and complements Article I, in which simulations and experiments on gA in mixtures of short and long lipids are reported. The results indicate that the concentration of lipids that match the length of the dimer is enriched in the first shell. From the experiment, the complete dimer has a longer lifetime in the mixture, indicating that it is stabilized by the matching lipids. From the simulation, the concentration of the lipids can be directly observed. The results of the modeling described in this article both provide an estimate of the energetics of lipids that drive redistribution and present independent evidence for the behavior of the lipid-channel boundary condition itself. The results are in excellent agreement with the computed and measured redistribution from Article I. Results taken from Article I are as follows: 1) Lipid redistribution timescales, to build a kinetic model of lipid redistribution necessary to estimate simulation error bars; 2) Experimental and simulation lipid redistribution, to provide independent validation of the boundary condition determined by the model; 3) Lipid shell counts, to estimate the compression energy of long lipids near the short channel, This study provides theoretical support for Article I in the following ways. 1) It gives estimates of the compression energy of the long lipids into the first lipid solvation shell of gA. 2) It gives estimates of the standard deviation of the enrichment expected given the timescales from Article I and of the energetics in this work. 3) It holds implications for the structure of the lipids near the channel, given the observed redistribution. The general effectiveness of this model beyond gA will be determined by the limit of its applicability; the first place to look is in the first shell of lipids around the protein, where the lipid material ends and the complexity of the protein begins. As part of the introductory material for this work, the theoretical background for first-shell lipid properties, including compression and curvature stress, will be discussed. Table 1 lists descriptions of the symbols and abbreviations used frequently in this work. Leaflet bending and compression For a system with an inclusion (like gA) whose hydrophobic thickness does not match the surrounding lipid environment, the soft leaflet must deform to match the length of the inclusion. This article deals with the case where leaflets are

1200 Biophysical Journal 112, 1198–1213, March 28, 2017

TABLE 1

Summary of Symbols and Abbreviations

EL kc;m KA;m K A ðzÞ H d 0

F ð0Þ

c0 cr ct hðrÞ r0 h0 hp hns t0 2D-CEM(s ¼ 0) 2D-CEM(s free) 3D-CEM

Surface model lipid deformation energy (kcal mol1) Single-leaflet (monolayer) bending modulus (kcal mol1) Single-leaflet area compressibility modulus ˚ 2) (kcal mol1 A Local compressibility modulus per unit thickness ˚ 3) (kcal mol1 A Phenomenological model force coefficient for gA ˚ 2) hydrophobic matching (kcal mol1 A Difference in hydrophobic thickness between the ˚) channel and the TS for dissociation (A Simulation measure of curvature stress: the derivative of the free energy with respect to curvature, evaluated ˚ 1) at zero curvature (per unit area, kcal mol1 A 1 ˚ Single-leaflet spontaneous curvature (A ) ˚ 1) Radial value of the local principal curvature (A Tangential value of the local principal ˚ 1) curvature (A The hydrophobic thickness of the leaflet surface in ˚) the 2D-CEM models (A ˚) Effective radius of the gA channel (A Average hydrophobic thickness of a lipid leaflet with ˚) no peptide (A Effective hydrophobic thickness of the gA ˚) monomer (A Position of the neutral surface of bending of the ˚) leaflet (A Average total thickness of a lipid leaflet with no ˚) peptide (A Surface continuum model (two-dimensional), with boundary condition fixed to zero slope Surface continuum model (two-dimensional), with boundary condition free to minimize the leaflet energy Solid material continuum model (three-dimensional), with boundary condition fixed by all-atom simulated ‘‘shape’’

thicker than the inclusion, and so the language will follow this case; given the strong energetics of hydrophobic matching, lipids in the first solvation shell must be compressed in thickness to match the channel. Following the deformation outward radially, the leaflet must eventually bend up to match the thickness of the lipid environment far from the inclusion. However, the precise behavior of this deformation in the interim must be modeled or inferred from experiment (30,31,44). The widely used 2D-CEM (30,31,34,43–45) combines the Helfrich Hamiltonian (4,5,7) with the modulus for lipid area compressibility to account for curvature (here, in the Monge gauge (46)) and compression, respectively, through the function describing the leaflet hydrophobic surface, hðzÞ, " Z 2 kc;m
EL ¼ dA  V2 hðrÞ  c0 2 2 #  KA;m hðrÞ  h0 ; (2) þ 2 h0

Modeling the Gramicidin-Lipid Boundary

where EL is the deformation energy, kc;m is the bending modulus (8), KA;m is the area compressibility (3,47), and c0 is the spontaneous curvature (19) of the lipid matrix. Experimentally, KA;b ðKA;b ¼ 2KA;m Þ can be determined by the pipette aspiration technique of Evans to be between 265 mN/m at 291 K (47) and 290 mN/m (48), depending somewhat on the application of tension. The work of (47) relates kc;b to KA;b through the ‘‘polymer brush’’ model as kc;b ¼

KA;b ð2h0 Þ2 ðpolymer brushÞ; 24

(3)

where the subscript b indicates a bilayer quantity. Pan et al. (8) test this model using x-ray scattering to determine that kc;m ¼ 6.2 kcal/mol and h0 ¼ 27.7 at 288 K for DOPC. From their kc;m, the inferred value of KA;b (260 mN/m) matches the pipette aspiration value. An isotropic bilayer continuum model with uncoupled leaflets yields kc;b ¼ ðKA;b ð2h0 Þ2 =48Þ. The assumption in this work is that the acyl chain region of the bilayer is mechanically softer, accounting for the difference in the denominator. This is equivalent to fitting the local elastic moduli of the chains to recover the polymer-brush relation. An all-atom simulation test of Eq. 3 using, to our knowledge, new methodology to determine kc;m found good agreement for DOPC (KA;b ¼ 290 5 20 mN/m, kc;b ¼ 17 kcal/mol) and similar lipids, although fully saturated chains were outliers (40). If the material is chosen to be inhomogeneous, as in this work, the local moduli can be fit to the observed experimental scaling. These data are consistent for KA;m and within 25% for kc;m. However, there are a variety of techniques for measuring kc;b , including fluctuation analysis of vesicles (49). Fluctuation analysis in particular shows a wide range of kc;b (50), perhaps as a result of buffer conditions (51). This work uses CEMs with kc;m and KA;m matched to the simulation values, for comparison with curvature stress and redistribution simulations. Previous modeling has pointed out that the total energetics of the deformation depend critically on the behavior of the deformation in the first shell, frequently described using the slope of the deformation as it meets the inclusion (31,52–54). In the Appendix, the energetics of compression and curvature are analyzed with relation to the boundary condition. The well-known relation between the boundary condition and total curvature stress is repeated, as is an estimate of how the slope influences first-shell compression. Two points pertain to the validation metrics discussed in the Introduction. First, the total curvature stress of the leaflet is proportional to the value of the slope. Second, the compression of lipids in the first shell is very large for small slopes, where the leaflet is not allowed to immediately deflect away from the inclusion in a way that would minimize compression. This first shell compression will lead to a large energetic penalty for long lipids to occupy the first shell. These effects are illustrated in Fig. 1. Through these two

FIGURE 1 The consequences of the slope, h0 ðr0 Þ, of the membrane surface intersecting the inclusion, that relate the first-shell compression and curvature stress. If the slope is constrained to be zero, there is zero net curvature but large compression (red). If the slope is unconstrained, there is net positive curvature (blue). See Eqs. 36 and 37. To see this figure in color, go online.

metrics, the simulations and analysis below justify a model that considers the 3D deformation of the leaflet around the inclusion, rather than simply a 2D model of the hydrophobic surface. An important benefit of the 3D-CEM, with information from molecular simulation, is the removal of the ambiguity of the boundary condition. As shown below, although not quantitatively as accurate for curvature stress, the properties of the 2D-CEM(s free) are comparable to those of the 3D-CEM. However, without molecular justification for the free boundary condition, extension to larger proteins is questionable. The blue and red arrows in Fig. 1 illustrate the slope boundary condition near the peptide for positive and zero slopes, respectively. In this figure, the zero slope in red shows the correlation with increased compression. The relation between slope, compression, and curvature is summarized in the figure by color. These experimental and simulation observables (curvature stress and compression) are critical to assessing model quality. Three models are presented below. The first two are continuum models that reduce the effect of the bilayer to that of a surrounding material. These models are useful for rapidly predicting the effect of the surrounding lipids on the inclusion. The last model is the all-atom model, which derives its predictive power from the molecular detail. For the all-atom model, the deformation energy itself is ambiguous, yet the curvature stress can be directly calculated. As discussed above, curvature stress is directly related to the amount of first-shell compression. The Results section is structured to first validate the curvature stress of the 3DCEM against that of the all-atom model, and then to use

Biophysical Journal 112, 1198–1213, March 28, 2017 1201

Sodt et al.

the 3D-CEM to predict the degree of lipid redistribution around the channel, which is reported in Article I. MATERIALS AND METHODS 2D-CEM A 2D-CEM is defined by considering Eq. 2 directly. In this model, the deformation of the bilayer normal to the surface has been integrated implicitly, that is, the local curvature and height at any point on the surface defines the local lateral deformation away from the surface. For the case of a cylindrical inclusion, the lack of flexibility for the off-surface deformation leads to only one choice for h0 ðr0 Þ that is consistent with the shape of the protein, h0 ðr0 Þ ¼ 0. Within this framework, deviations from the ideal cylindrical shape may lead to non-zero values of h0 ðr0 Þ. The choice of h0 ðr0 Þ to match simulation profiles is discussed below. However, two models are chosen as bookends to encompass the observations of the simulations. The first, denoted 2D-CEM(s free), lets the slope of the membrane at the inclusion boundary vary to minimize the deformation energy. The second, denoted 2D-CEM(s ¼ 0), constrains the slope of the membrane at the inclusion boundary to be zero, consistent with a cylindrical inclusion. A more complex energetic description of the slope near the inclusion is possible, for example, by fixing the value of the slope to match a molecular simulation. This has proved to be difficult, however, as shown in (54), where, even for a thick bilayer around gA, the upper leaflet appears to slope down before sloping back up. The use of the hydrophobic surface to extract the slope leads to negative values of h0 ðr0 Þ, leading to negative curvature through Eq. 38. Additionally, this reversed slope would lead to a very strong contribution to the deformation energy in the first shell; in Eq. 41, ffs would be negative for a thick bilayer.

2D-CEM methods The energy of the 2D-CEM is computed using a piecewise C1 continuous cubic spline curve representation of hðrÞ. Twelve equally spaced radial control points fix the value of the spline, with continuity of the first derivative acting as a restraint. For the inner boundary condition, the value of hðrÞ is fixed to the inclusion thickness. The slope of the membrane as it meets the peptide can be fixed or optimized, determining the two models (s ¼ 0 or s free). An additional restraint is placed on the lipid density. The hydrophobic volume beneath hðzÞ is computed and restrained with a very high force constant to its undeformed value such that it effectively acts as a constraint. The restriction to a finite area (equivalent to that of the simulation) and the use of a spline result in small differences compared to analytical approaches. The values of the control points are then optimized via the downhill simplex (55) method until the energy (Eq. 2 plus the volume restraint) is minimized. In both the 2D-CEM and the 3D-CEM (described below), the radius ˚ and the single-leaflet thickness is set to 13 A ˚ . Both of of gA is set to 10 A these quantities are ambiguous; as they are resolved at the atomic level, they must be set phenomenologically. In this work, they are set to best represent ˚ for the singlethe deformations seen in Fig. 3. Note that our use of 13 A ˚ (56) inferred leaflet thickness of gA differs from the value near 11 A from the small change in peak-to-peak distance in x-ray scattering of gramicidin in dilauroylphosphatidylcholine at a 1:10 ratio of peptide to lipid, versus pure dilauroylphosphatidylcholine.

values of KA;m and kc;m defined for the surface to local moduli of the entire thick material. This allows a consistent treatment of the origins of bending and compression, including their coupling away from the surface. This naturally happens at the expense of introducing more parameters, as described below. From the perspective of the bilayer as a 3D material with local elastic properties, the area compressibility modulus of a leaflet can thus be viewed as the average of the local material area compressibility modulus:

Z

KA;m ¼

0

t0

dzK A ðzÞ;

(4)

where t0 is the thickness of the leaflet and the units of K A ðzÞ are those of pressure (that of an elastic modulus of a material in three dimensions). A bending deformation induces a lateral deformation proportional to the distance from the neutral surface (e.g., the area of a cylinder is proportional to its radius). The local lateral deformation energy depends on the square of the deformation through K A ðzÞ. The bending modulus is thus:

Z

kc;m ¼

0

t0

dzK A ðzÞðz  hns Þ2 ;

(5)

where the z position of the neutral surface (58) is hns . It is defined by

Z 0

t0

dzðz  hns ÞK A ðzÞ ¼ 0:

(6)

In the 3D-CEM, each tetrahedral volume element of the leaflet is assigned a local area compressibility modulus, K A ðzÞ, according to its distance from the bottom of the leaflet. The quantities KA;m , kc;m , and hns are experimentally available and thus can parameterize at most three parameters of the function K A ðzÞ through Eqs. 4–6. A simple choice is three regions, each with a different constant value of K A ðzÞ. The thickness of the regions is defined by the varying chemical elements: lipid acyl chain, glycerol/ headgroup, and hydrated headgroup. With varying acyl chain length, the thickness of the acyl chain region is extended without changing the value of K A ðzÞ. The value of K A ðzÞ defines local elastic moduli for the volume element:

K A ðzÞ ¼ lxx;xx  2lxx;zz þ lzz;zz ;

(7)

where lxx;xx ¼ lzz;zz is the bulk modulus (18). The bulk modulus used is ˚ –3, consistent with the bulk modulus of alkane material 0.125 kcal/mol/A (59), and which is much larger than K A ðzÞ, thereby allowing only small fluctuations around the average density (45). In addition, the modulus ˚ –3 controls the penalty to lateral shear, a paramlxz;xz ¼ 0:00537 kcal/mol/A eter that is set according to the phenomenological lipid tilt modulus (60). The dimensions x and y are equivalent in the continuum model, and thus, all moduli are unchanged by exchanging x with y. This is consistent with lateral fluidity. The values of K A ðzÞ are tabulated in Table 2. To illustrate the connection between the moduli and the deformation, an example will be considered. The deformation field uðrÞ is optimized to parameterize deformed coordinates r0 that are compatible with the boundary conditions:

r0 ¼ r þ uðrÞ:

(8)

If the undeformed width of an element is w, that is, the x values of the left (l) side and right (r) side of the box are xl and xl þ w:

3D-CEM In Eq. 2, KA;m and kc;m are available from experiments. However, these experiments are only sensitive to large-lengthscale deformations of the bilayer, giving rise to the ambiguity of the deformation in the first shell. The 3D-CEM, developed for the membrane by Campelo et al. (57) to predict the curvature induction of membrane-embedded peptides, relates the

1202 Biophysical Journal 112, 1198–1213, March 28, 2017

w ¼ xr  xl ðundeformedÞ;

(9)

the deformed width will be

w0 ¼ xr  xl þ ux ðxr Þ  ux ðxl Þ ¼ w þ ux ðxl þ wÞ  ux ðxl Þ (10)

Modeling the Gramicidin-Lipid Boundary TABLE 2 Values of the Local Area Compressibility Modulus Used in the 3D-CEM K A ðzÞ ˚) Lower Limit (A 0 h0 h0 þ 7

˚) Upper Limit (A

˚ 3 kcal mol1 A

mN m1 nm1

h0 h0 þ 7 h0 þ 11

4.28  103 20.0  103 8.25  103

35.6 97.2 23.0

ture at the hydrophobic surface would imply that the area is increased near the tails. However, this may be incompatible with the shape of the inclusion. In this work, the boundary condition at the lipid-gA interface is extracted from the molecular dynamics simulation. Shown in purple in Fig. 3 are the hydrophobic surfaces of the leaflet for gA þ dC18:1 (top) and gA þ dC22:1 (bottom) simulations. Also shown are traces of the path of lipids, which effectively indicates the local director of the leaflet. Local lipid directors are calculated as reported in the Materials and Methods (allatom) section.

Lipids of varying length are modeled by changing h0 .

3D-CEM methods

and

w0  w ux ðxl þ wÞ  ux ðxl Þ ¼ : w w

(11)

The 3D-CEM energy is computed as a sum of energies of each local finite element of the membrane. The energy of a finite-volume element is computed as (18)

Z

In the limit of a small volume element, w, x ðxl Þ lim ux ðxl þwÞu w w/0

¼ uxx ðxl Þ;

(12)

where uxx ðxl Þ is the derivative of the x component of uðxÞ with respect to x evaluated at xl . Thus,

w0  w w0 ¼  1zuxx ðxl Þ: w w

Ei ¼

(13)

V

1 dr3 lij;kl uij ðrÞukl ðrÞ; 2

(15)

where uij ðrÞ is the jth derivative of the ith Cartesian component of the deformation gradient field, uðrÞ. Each vertex a of a tetrahedron is an optimized variable of the model. Neighboring tetrahedra share vertices to maintain material continuity. The value of uðrÞ is defined at each of the vertices of the tetrahedron by subtracting its position from its initial reference position:

As discussed further below, the elastic model for this deformation energy density is

ua ðra Þ ¼ ra  r0a :

1 2 E ¼ lxx;xx uxx ðxl Þ : 2

The field uðrÞ is then linearly interpolated using a barycentric coordinate system to evaluate the tetrahedron energy using Eq. 15. The set of vertices that neighbor the modeled inclusion are constrained to adhere to the shape of the inclusion. Points at the bilayer midplane are constrained to be at z ¼ 0, modeling leaflet symmetry. Additionally, the vertex defining the border between the acyl chains and headgroups is constrained to be at the channel hydrophobic thickness, modeling the strong condition of hydrophobic matching. The inclusion is modeled as an ‘‘inner’’ boundary condition, where vertex points are constrained to lie on the surface, parameterized by a spline curve, rðzÞ, that approximately fits the local director as extracted from the all-atom simulation. The curve is defined as

(14)

This element is illustrated in Fig. 2, along with illustrations of gA in the deformed model alongside the all-atom model. If constant volume is (approximately) maintained, uzz ðrÞz  uxx ðrÞ, and Eq. 7 is recovered. The 3D-CEM is able to model deformations away from the hydrophobic surface that are consistent with the molecular shape of the inclusion, previously shown to be important for bilayer deformations (36,61). For example, in the 2D-CEM, the function hðrÞ does not provide details for how the material changes near the ends of the acyl chains. Rather, negative curva-

rðzÞ ¼

8 > > > > > < > > > > > :

r0 þ rc

r0 "qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

1  ðz  h0  DÞ2  1

r0 þ rc ½sinðq0 Þ  1 

z  zq tanðq0 Þ

#

(16)

if z%h0  D if h0  D < z < h0  D þ zq

;

(17)

if z > h0  D þ zq

FIGURE 2 Illustration of the basics of the 3DCEM using the system gA þ dC24:1 as an example. The relation between the width, w, of a volume element and energetics is discussed in relation to Eq. 14 in the main text. To see this figure in color, go online.

Biophysical Journal 112, 1198–1213, March 28, 2017 1203

Sodt et al. of the simulation. The lipids diffuse on the surface as it fluctuates according to the simulated ensemble, and the deformation can be computed according to the average value of a particular site in the model chosen to correspond to the hydrophobic surface. In theory, the molecular model contains all of the relevant energetics of the surface. However, they cannot be extracted trivially. In particular, the energetics are not a local function of the height or material deformation, as they are in the continuum models, and so a local energy density as in the integrand of Eq. 2 cannot be computed. Even with the quantity hðrÞ computed from a histogram of atomic positions, the inclusion-lipid boundary condition is still ambiguous, as can be seen from the hðzÞ values plotted in purple in Fig. 3. The quantity hðrÞ is defined even above the peptide, where the surface is not behaving as a continuum membrane. These must represent fluctuations of the hydrophobic atom above the peptide. For this case, Eq. 2 cannot be applied to hðrÞ in the vicinity of the peptide. Interpreted directly, they indicate that depending on where the slope is radially evaluated, it may be zero or even negative (54). The trace of the lipid atom positions, representing the local director of the lipid traced along its molecular path, is an alternative representation of the boundary condition and shows a dramatically different value of the slope (i.e., the slope at a point in z is defined as a radial vector orthogonal to the trace) that is positive for reasonable choices of the hydrophobic surface. In Article I, extensive (multiple microseconds) simulations of gA in mixed lipid bilayers are reported. In these simulations, short acyl-chain lipids are shown to reside preferentially in the first shell around inclusion, but with a fraction indicating only moderate lipid compression. In addition, the work reports the timescale of lipid solvation shell transitions. Both the occupancy and timescales are used later in this article to assess model energetics.

All-atom model methods

FIGURE 3 Lipid deformations from the all-atom molecular model. The purple line represents the hydrophobic surface computed from the average position of the acyl chain carbon near the lipid surface. Black lines represent traces of lipids, computed as described in the main text. To see this figure in color, go online. where rc , D, and q0 parameterize the surface, and zq ¼ rc cosðq0 Þ. The function rðzÞ is constant up to the inclusion thickness (minus D), then curves inward with radius rc through angle q0 , then proceeds linearly with continuous slope. The deformation around the inclusion may be quite large. Thus, the ‘‘Green-Lagrangian’’ finite-strain tensor is appropriate. With this measure of strain, the coefficients uab are modified as

u0ab

! X vug vug 1 vua vub þ  : ¼ 2 vb va g ¼ fx;y;zg va vb

(18)

Additionally, in the vicinity of the peptide where the strain is the largest (within ˚ of the peptide boundary at r0 ), the field uðrÞ is required to be cylindrically 8A symmetric. Together these constraints, including the use of periodically replicated simulation cells, define the boundary conditions of the model under which the energy is minimized using a standard BFGS (62) algorithm.

All-atom model The molecular model of the leaflet deformation does not explicitly account for curvature and compression, rather, they are present naturally as a result

1204 Biophysical Journal 112, 1198–1213, March 28, 2017

The initial coordinates for all systems were obtained from the Membrane Builder module in CHARMM-GUI (www.charmm-gui.org) (63–65), and the gA structure was from PDB: 1JNO (66). Pure bilayer simulations were conducted of phosphatidylcholine (PC) lipids with varying tail length, dC16:1, dC18:1, dC20:1, dC22:1, and dC24:1, with 50 lipids/leaflet, 4550 waters/lipid, and 0.15 M KCl. Two sets of gAþbilayer systems were made for each lipid type, a set with two gA monomers and a set with a single gA dimer. Systems containing gA were built with 90 lipids/leaflet, ~45 waters/lipid, and 0.15 M KCl. To eliminate differences in gA structure due to fluctuations, gA in all monomer and dimer systems were restrained by a harmonic root mean˚ 2) square deviation (RMSD) restraint with a force constant of 5.0 kcal/(mol$A relative to the minimized PDB: 1JNO structure. In all cases, the gA was parameterized using the C22 protein parameter set with the ‘‘correctionmap’’ (CMAP) addition (67), and lipids were parameterized using the C36 parameter set (39). Bonds with hydrogen atoms were constrained using the SHAKE algorithm (68). To reduce possibly long equilibration times, bilayers of systems including a gA dimer were deformed to an idealized surface from a Helfrich continuum model. The surface itself was determined by a deformation taken from the approximated length and radius of gA. By doing this, lipid headgroups could be placed close to the hypothetical equilibrium conformation before any minimization was performed. Monomer gA systems contained one monomer per leaflet at maximal separation (to minimize transmembrane monomer-monomer interactions) and with the z position taken from previous results (53). Monomer x and y positions were held by ˚ ), and were a harmonic potential with a force constant of 1.1 kcal/(mol$A unconstrained in z. No positional restraints were placed on the dimer. All simulations (pure, dimerþbilayer, and monomerþbilayer) were run with a 1 fs time step for 120 ns using NAMD. Constant pressure of 1 bar was regulated by Nos
e-Hoover Langevin piston pressure control (piston period of 50.0 fs and piston decay of 25 fs), and constant temperature of 310.15 K was maintained by Langevin dynamics (coupling coefficient

Modeling the Gramicidin-Lipid Boundary of 1 ps–1). Electrostatic interactions were calculated using the particle-mesh ˚ ; k ¼ 0.34 A ˚ –1; and sixth-order B-spline Ewald method (69) (mesh size, 1 A interpolation), van der Waals interactions were switched off between 10 and ˚ by a force-switching function (70), and a pair list distance cutoff of 12 A ˚ was imposed. All but the PME contribution to the pressure profile 16 A was calculated during simulation using NAMD’s pressure profile capability (71,72), with 100 equally spaced slabs and the profile calculated every 100 fs. The PME contribution was calculated post-simulation, where the contribution was calculated every 0.5 ps.

radius R1 ¼ cglobal ¼ 0. Consider a small deformation of the membrane to adapt to the inclusion boundary condition that changes the curvature at a particular point to ec . This deformation can be ‘‘undone’’ by changing the global curvature to ec , as global curvature is position independent (a cylinder), and small curvature deformations are additive. The system 0 will appear to have spontaneous curvature ec , so from Eq. 22, F ð0Þ is evaluated as

Computing curvature stress

As a check on this approximation, the numerical procedure for evaluating the 2D-CEM energy was modified by introducing two transformations to the coordinates for evaluation: first, the ‘‘tension’’ transformation is applied,

The curvature stress, represented as the derivative of the free energy per unit area ðFÞ with respect to curvature, is computed in all three model types, the 2D-CEM, the 3D-CEM, and the all-atom, using the methodology defined  0 0 below. In each case, it is computed as F ðR1 ¼ cÞ  c¼0 ¼ F ð0Þ, where R1 ¼ c is the global curvature of the system (e.g., the system is embedded in a cylinder of radius R).

0

The value of F ð0Þ can be computed from a 3D model (including the allatom) using the first moment of the lateral pressure profile (17,71–74):

Z

0

dz z½pT ðzÞ  pN ðzÞ:

(19)

0

Equivalently, F ð0Þ can be computed by the finite difference: 0

Z

dr2pr½cr ðrÞ þ ct ðrÞ

ut ðrÞ ¼ fexpðet Þx; expðet Þy; expð2et Þzg;

(24)

(25)

where et is an independent variable that is optimized to minimize the energy and assures zero tension. Then, the finite-difference derivative with respect to ec , the global system curvature relative to the z axis, can be taken:

  1
uc ðrÞ ¼ ec xz; yz;  x2 þ y2 þ 2z2 : 2

0

F ð0Þ computed for the 3D-CEM and all-atom models

F ð0Þ3D ¼ 

0

F ð0Þ2D ¼ kc;m A1

(26)

A few modifications that complicate the algorithm were necessary. Leaflet thickness was computed relative to the bilayer midplane that was modified by the same transformation. Hydrophobic mismatch and protein compressibility were changed from constraints to nearly equivalent stiff harmonic restraints, but which could be modified by the transformations. The modi0 fied algorithm produces values of F ð0Þ equal to Eq. 24 for the Monge gauge.

1

F ð0Þ ¼ ð2eÞ hVðfr þ ue ðrÞgÞ  Vðfr  ue ðrÞgÞi; (20) where



 1
2 2 ue ðrÞ ¼ e xz; 0;  x þ z : 2

(21)

Here, ue ðrÞ preserves density, and so the kinetic energy component of the free energy cancels in molecular models, with an important caveat regarding constraints (e.g., SHAKE) discussed in (71). The all-atom model has the added complexity that an individual monomer of gA deforms the leaflets and so induces curvature stress even when it does not form a transmembrane channel. This is accounted for by performing simulations of both the monomer and dimer systems and 0 computing F ð0Þ ¼ Fdimer  Fmonomer under the same peptide/lipid ratio, where the subscript indicates the molecular simulation performed. This expression attempts to account for the change in curvature stress due only to the leaflet deformation, the closest comparison possible to the continuum models.

Evaluating shell energetics as a function of lipid length To estimate the effect of lipid redistribution from a CEM, it is necessary to estimate the effect of replacing a lipid of a certain length with another lipid species with a different hydrophobic length. This is accomplished by obtaining a model hðrÞ (2D) or uðrÞ (3D) from a model using the average lipid length of a mixture, then modeling what would happen if a lipid with different thickness were placed at r. For the 2D-CEM, the approach is straightforward: the energy per lipid for an alternate lipid with height h00 is

Eshell ¼

A0 2 pðr1  r02 Þ

Z

r1 r0

0

Curvature in the 2D-CEM is accounted for locally by the Helfrich Hamiltonian,

FðcÞ ¼ 0

kc;m ðc  c0 Þ2 2

F ð0Þ ¼ kc;m c0 ;

(22) (23)

where F is the local energy per unit area at a point with curvature c. The global curvature stress is approximated by computing the integral of the local curvature of the system that is embedded in a cylinder with inverse

2  KA;m hðrÞ  h00 dr: 2 h00

(27)

In the 3D-CEM model, information beyond hðrÞ is available, namely, the length of the deformed lipid taking into account local curvature:

Z heff ðrÞ ¼

F ð0Þ computed for the 2D-CEM

2pr

h0 0

dzð1 þ uzz ðrÞÞ

(28)

where heff ðrÞ is the effective height of the lipid at the hydrophobic surface taking into account the deformed path. The integral is performed over a path through z defined in the undeformed state, so h0 is the average height of the mixture. The quantity heff ðrÞ is then used in lieu of hðrÞ in Eq. 27. The radii r1 and r0 in Eq. 27 are chosen to select a particular shell of lipids. These values are taken from Article I. Additionally, the compression energy of the lipid in the bulk, beyond the fourth shell, is used as a reference and is thus subtracted from these values. The effect of curvature is ignored; it is assumed that curvature does not change significantly as a function of lipid redistribution and so the energy is invariant. Even if curvature does change, the majority of the energy is from compression.

Biophysical Journal 112, 1198–1213, March 28, 2017 1205

Sodt et al.

Local lipid director (trace) These traces are constructed by first accumulating a histogram (on a fine mesh) of both the density, ri ðrÞ, and average height, zi ðrÞ, of each atom site i of the lipid as a function of the r, the distance from the peptide. Each trace averages over n lipids. Second, the value r0 is computed such that

Z

r0 0

drri ðrÞ ¼ n;

(29)

that is, such that n lipids are accounted for between 0 and r0 . Then, the average radius and height of the lipid site are computed as

fhri i; hzi ig ¼ n

1

Z

r0

0

drri ðrÞfr; zi ðrÞg:

(30)

These sets of fhri i; hzi ig pairs are then connected along the chemical connectivity of the lipid to yield the trace. The vector between sites fhri i; hzi ig and fhriþ1 i; hziþ1 ig then defines a locally averaged lipid director. For the next trace, r0 is used as the lower limit, and the upper limit of integration is computed to account for the next n lipids. Continuous connectivity is achieved by using only the sn-2 chain and ignoring atoms that branch off the path to the headgroup.

Monte Carlo modeling of lipid diffusion For each selected frame of an all-atom simulation, lipids were assigned an average position and a 2D Voronoi decomposition was computed. Sites that share a border were labeled as ‘‘connected,’’ establishing a map of neighboring sites. Each site was assigned a shell based on the minimal number of neighbors to cross to reach the inclusion. The total energy for the model runs is computed by assigning each of the short or long lipids (labeled i) an energy ðej Þ by shell j, where sij ¼ 1 if lipid i is in shell j, else sij ¼ 0:

E ¼

X

ej sij :

(31)

i

The shell energies are extracted from either the 2D- or 3D-CEM. For the Monte Carlo simulations, the initial lipid map was used to establish neighbors for the ensuing model run. Lipid hops for each lipid are attempted every 0.2 ns, a period much shorter than the lipid hop time (t ¼ ~70 ns). The probability used was kð0:2=tÞ. The value of k was adjusted to be 0.8 to match the simulation value of t. For this attempt, the lipid and one of its randomly selected neighbors switch places, and the energy is recomputed and accepted or rejected with the standard Metropolis criterion (75). For both the Monte Carlo simulation and the values extracted directly from the simulation, the enrichment of the populations of each shell for each long or short lipid chain is computed.

RESULTS AND DISCUSSION Deformation The hydrophobic surfaces, hðrÞ, of the four models are shown in Fig. 4 for gAþCd22:1 (top) and gAþCd24:1 (bottom). Plots for gAþCd16:1, gAþCd18:1, and gAþCd20:1 are shown in the Supporting Material. For the all-atom model, the CHARMM ‘‘C21’’ atom is plotted (the ester atom on the sn-2 chain). For the CEMs, hðrÞ is plotted for the interface between the polar and apolar regions where the matching condition is chosen. The key feature of the plot is that

1206 Biophysical Journal 112, 1198–1213, March 28, 2017

FIGURE 4 Hydrophobic surfaces plotted for the all-atom model and CEMs. The plots are color-coded as follows: purple for the all-atom model, red for the 2D-CEM(s ¼ 0), blue for the 2D-CEM(s free), and green for the 3D-CEM. To see this figure in color, go online.

the 2D-CEM with slope unfixed and the 3D-CEM match the all-atom simulation reasonably well, whereas the 2DCEM with the slope fixed at zero is qualitatively incorrect. 0

Comparison of F ð0Þ between models 0

As noted by Eq. 23, F ð0Þ is proportional to the total curvature of the system and so should be sensitive to the effect of the boundary condition (directly through Eq. 38 for the 2D0 CEM). Results for F ð0Þ computed for all models are shown in Fig. 5. For both 2D-CEMs, the total curvature (in the Monge gauge) is equal to Ring’s expression, 2pr0 s. Thus, 0 the s ¼ 0 model has F ð0Þ ¼ 0 for all lipids, inconsistent with the simulation results. With s free, the value is qualitatively similar in that the all-atom simulation and 2D elastic model both indicate net positive curvature. The 3D-CEM values are consistent with the all-atom simulation. In two ways, the 3D-CEM is a more faithful model of the allatom simulation. First, the 3D-CEM boundary condition has been chosen to match the all-atom simulation. As discussed, in the 2D models, curvature stress is directly dependent on this boundary condition, and the same is true of the 0 3D-CEM. Second, calculation of F ð0Þ from the lateral stress profile in the 3D-CEM is equivalent to calculation in the all-atom model. Even without interpretation through 0 the continuum models, the values of F ð0Þ from molecular

Modeling the Gramicidin-Lipid Boundary

0

FIGURE 5 Global curvature single-leaflet stress, expressed as F ð0Þ. Points are evaluated from the all-atom simulation. For the 2D-CEM, thin lines use the Monge gauge to evaluate curvature, whereas thick lines use the true definition. To see this figure in color, go online.

simulation suggest that the s ¼ 0 boundary condition model is incorrect even without being able to unambiguously characterize the slope of the all-atom model hydrophobic surface. That is, even though a strict evaluation of s from the experimentally deduced value of H may suggest a particular value of s, this is likely to be problematic and to reflect the general difficulty of applying continuum models to describe molecular events. The peptide-lipid boundary In the pioneering study of Huang (30), the shape of the bilayer deformation around a gramicidin channel in monoglyceride membranes was deduced from the variation in single-channel lifetimes as a function of bilayer thickness. Using the available elastic moduli for monoglyceride bilayers, the experimental results could be accounted for by a 2D-CEM only if s was negative, quite different from the positive value of s that would minimize the deformation energy (29,31). Invoking the minimum-energy principle, this could imply that the system is not at equilibrium, which is difficult to reconcile with timescales of molecular equilibration,  100 ns, found in Article I. Alternatively, one or more of the assumptions involved in evaluating H, and thus s, may be incorrect. The elastic moduli close to the channel may, for example, differ from the continuum moduli that are used in the 2D-CEM (76), though this is difficult to reconcile with the more recent molecular dynamics simulations of Kim et al. (53). There could be unanticipated energetic contributions to the deformation energy that are not included in the 2D-CEM, such as from tilt (31) or the Gaussian curvature modulus (22,45). However, tilt is incorporated into the 3D-CEM, and the latest estimates of the value (10,60,61) do not sufficiently raise H

for the s free model to impact the results of Huang. Gaussian curvature is discussed below. The difficulty linking H and s may simply reflect the well-known limitations associated with the strict application of continuum models to describe molecular events, which in turn implies that H should be regarded as a phenomenological parameter describing a particular experimental system. This study is concerned with bilayers formed from phospholipids with two acyl chains, and so a direct correspondence with monoglycerides cannot be expected. Nevertheless, the values of H extracted from lifetime measurements in monoglycerides and phospholipids are similar (69 kJ/mol/nm2 (38) and 56 kJ/mol/nm2 (34), respectively), although the phospholipid system is more challenging, requiring added decane. Bilayers formed from phospholipids, however, have bending and area compressibility moduli that are two to three times higher than those formed from monoglycerides. H for gramicidin channels in phospholipid bilayers thus would be expected to be correspondingly higher, given similar boundary conditions—though the possibility cannot be excluded that the boundary condition may shift as a result of mechanical parameters or the molecular details of the lipids. A picture thus emerges of the difficulty in relating lifetime measurements to the first-shell behavior of lipids inferred from the continuum model. The extracted boundary condition depends on accurate values for the mechanical parameters, the effect of added solvent such as decane, and qualitative differences in the lipids themselves (phospholipids versus monoglycerides). Acknowledging these difficulties, the 2D-CEM(s free) and 3D-CEM values of H (Table 3) are in reasonable agreement with previous estimates of 56 kJ/mol/nm2 for hydrocarbon-containing bilayers (34). Both sets of values are less than the 2D-CEM (s ¼ 0) value, 215.6 kJ/mol/nm2, but we consider the close similarity between the predictions of the 3D-CEM and the all-atom model (Fig. 5) to provide strong support for the 3D-CEM, subject to the limitations associated with applying continuum models to describe molecular events. Two, to our knowledge, new simulation techniques have been applied to extract the peptide-lipid boundary. First, the lipid-trace method revealed the form of the boundary lipids (Fig. 3) extracted from the molecular simulation. This boundary was used by the 3D-CEM to compute energetics and curvature stress. Second, the leaflet curvature stress was computed by all-atom simulation by subtracting 0 0 F ð0Þmonomer from F ð0Þdimer . The curvature stress of the TABLE 3 Phenomenological Deformation Energy Coefficients, H ˚ 2, kcal mol1 A ˚ 1 H kJ mol1 A

Model 2D-CEM(s ¼ 0) 2D-CEM(s free) 3D-CEM

215.6, 0.52 44.7, 0.11 35.7, 0.09

Biophysical Journal 112, 1198–1213, March 28, 2017 1207

Sodt et al.

3D-CEM matches well with the all-atom simulation, indi0 cating the importance of the boundary for implying F ð0Þ. The value of the phenomenological bilayer hydrophobic matching coefficient, H (defined in Eq. 1), is extracted from the models by computing the deformation energy under a series of matching conditions. The single-leaflet energy is multiplied by 2 (to model the total bilayer energy) and fit against the total bilayer mismatch (twice the sum of the single-leaflet mismatch). The coefficients for the three models are reported in Table 3. As expected, the 2D-CEM with slope fixed at zero has the largest value of H, as curvature is unable to reduce the large energetic penalty to firstshell lipid compression. The 2D-CEM(s free) and the 3D-CEM have relatively similar energetics. For the 3DCEM the energy is likely lower due to the more detailed treatment of first-shell compression; the path of lipid compression follows the bend of the lipid rather than being projected directly onto the z axis, where it will appear to be more extreme. Comparing a breakdown of curvature and compression energetics from the 2D-CEM(s free) model (where they are separable) indicates that ~75% of the energy is from compression and 25% is from curvature. The 2D-CEM(s free) and 3D-CEM values of H are in reasonable agreement with previous estimates for hydrocarbon-containing bilayers, 56 kJ/mol/nm2 (34), although the 3D-CEM value is somewhat too low. Three sources of error are likely. First, the experimental values are for hydrocarboncontaining membranes, which may preclude strict comparison. Second, the first-shell behavior of the 3D-CEM relies on extraction from the molecular dynamics simulation, which is subject to systematic errors in its parameterization. Furthermore, the mechanical description of these lipids, itself a combination of the boundary condition, spatially resolved elastic moduli, and assumptions of the continuum model, can legitimately be viewed skeptically. Third, the value of H extracted from lifetime measurements relies on the form of the TS (77). Increasing the estimate of the TS form implied by d directly decreases the inferred value of H. Model estimate of local redistribution The extracted model energetics are then used to predict the effect of lipid redistribution in a system with a mixture of short and long lipids. Values of the occupancy after ~3.5 ms molecular simulation of gAþdC24:1/dC16:1 and gAþdC22:1/dC18:1 mixed systems are used to assess the energetics of the CEMs. The model technique is to use a simple Monte Carlo simulation of lipid ‘‘hopping’’ between solvation shells that is determined by CEM energetics. The procedure for computing the energetics of mixtures for the CEMs is described in Materials and Methods. The compression energies from the three continuum models for the dC22:1 and dC24:1 lipids (used in the simulations and experiments of Article I) are reported in Tables 4 and 5, respectively.

1208 Biophysical Journal 112, 1198–1213, March 28, 2017

TABLE 4 Estimated Compression Energy of dC22:1 Lipid near gA, in gADdC22:1/dC18:1 Model 2D-CEM(s ¼ 0) 2D-CEM(s free) 3D-CEM

First Shell

Second Shell

Third Shell

0.49 0.33 0.27

0.29 0.13 0.06

0.11 0.04 0.00

All values are given in kcal/mol.

Summarized in Tables 6 and 7 are the occupancies of the first three shells after 3.5 ms of simulation of the four models, for the gAþdC22:1/dC18:1 and gAþdC24:1/dC16:1 mixed systems, respectively. The initial conditions for each model were set by the all-atom model. Although the redistribution predicted by the 3D-CEM is closest to the all-atom simulation, the difference between this and the 2D-CEM(s free) is not significant enough to warrant a strong statement that the 3D-CEM energetics are more accurate. However, the 2D-CEM(s ¼ 0) predicts significantly more rearrangement than in the all-atom model simulations. A few choices were made that affect the value of the enrichment, as summarized here. First, molecular simulations of the hydrophobic surfaces of gA in gAþdC24:1/ dC16:1, gAþdC22:1/dC18:1, and gAþdC20:1, described in Article I, show that the deformation surface is similar for each; they all have similar curvature. Thus, it is compression energetics that influence lipid enrichment in these mixtures. Second, simulations of gAþdC24:1/dC16:1 indicate that the tail density of dC16:1 in the first shell, which is similar to bulk, is shifted upward in the outer shells. This means that the lipid is not ‘‘stretched’’ when it moves to the thick bulk and so does not experience an energetic penalty. Thus, only the energetics of the tall lipids are considered, leading to lower enrichment. Third, we used mechanical constants (KA and kc ) that are compatible with only a subset of methods, which includes the all-atom force field. These choices do affect the calculated enrichment from the 3D-CEM, and were made with prior knowledge of the calculated enrichments from simulation and experiment in Article I. Dependence on bending modulus and Gaussian curvature modulus As noted previously, the bending modulus has the highest degree of experimental uncertainty. For the frequently TABLE 5 Estimated Compression Energy of dC24:1 Lipid near gA, in gADdC24:1/dC16:1 Model 2D-CEM(s ¼ 0) 2D-CEM(s free) 3D-CEM

First Shell

Second Shell

Third shell

0.70 0.47 0.43

0.45 0.21 0.17

0.20 0.07 0.01

All values are given in kcal/mol.

Modeling the Gramicidin-Lipid Boundary 0

TABLE 6 Shell Occupancy of dC18:1 Lipids in a Mixed gADdC22:1/dC18:1 Simulation Extracted over the Last 3 ms Model 2D-CEM(s ¼ 0) 2D-CEM(s free) 3D-CEM All-atom Experiment

Shell 1 0.66 5 0.62 5 0.59 5 0.50 5

0.03 0.03 0.03 0.04

Shell 2 0.58 5 0.02 0.54 5 0.02 0.51 5 0.02 0.51 5 0.02 0.58 5 0.04

Shell 3 0.50 5 0.49 5 0.49 5 0.52 5

0.02 0.02 0.02 0.01

For the models and experiment, values are represented as the mean 5 SE. Initial conditions and timescales were the same as for the all-atom simulation.

studied palmitoyloleoylphosphatidylcholine, x-ray scattering and pipette aspiration are typically on the lower end (kc;m ¼ 6.2 kcal/mol (78) and kc;m ¼ 6.5 kcal/mol (47)), whereas fluctuation analysis shows values approximately twice as high (50), apparently depending on buffer conditions. The CEM values in this article are affected if a larger or smaller bending modulus is used. For the 2DCEM(s free), the value of H decreases by 11% or increases by 27% if a bending modulus half or twice as large, respectively, is used. For the 2D-CEM(s ¼ 0) model, the value of H decreases by 27.6% or increases by 42%, respectively. Curvature stress for the 2D-CEM(s free) model is more sensitive and is effectively proportional to kc;m , with F0 ð0Þ reduced by 43.5% or increased by 81% if kc is halved or doubled, respectively. However, this variation is irrelevant to the conclusion that curvature stress validates (in part) the mechanical representation of the simulation boundary condition. The monolayer Gaussian curvature modulus, kg;m , is even more difficult to measure experimentally, and its impact on gA channels has not been well-characterized (45) but is hypothesized to contribute energetically to the bilayer deformation by Z EG ¼ 2p drrkg;m cr ðrÞct ðrÞ: (32) From analysis of phase diagrams and the Helfrich/Canham energy of methylated ethanolamine phospholipids, the value of kg;m is estimated to be 0.84 kc;m (79) for DOPC. If the term is added to the 2D-CEM(s free), the value of H inTABLE 7 Shell Occupancy of dC16:1 Lipids in a Mixed gADdC24:1/dC16:1 Simulation Extracted over the Last 3 ms Model 2D-CEM(s ¼ 0) 2D-CEM(s free) 3D-CEM All-atom Experiment

Shell 1 0.72 5 0.68 5 0.65 5 0.66 5

0.03 0.03 0.03 0.02

Shell 2 0.63 5 0.02 0.56 5 0.02 0.55 5 0.02 0.55 5 0.01 0.66 5 0.02

Shell 3 0.52 5 0.50 5 0.48 5 0.49 5

0.02 0.02 0.02 0.01

For the models and experiment, values are represented as the mean 5 SE. Initial conditions and timescales were the same as for the all-atom simulation.

creases by 67%, and the value of F ð0Þ is nearly completely ˚ . Unfortusuppressed, down to 0.022 from 0.099 kcal/mol/A nately, there is currently no straightforward way of including the hypothesized effect of kg in the 3D-CEM. CONCLUSIONS Intuitively, the information about the boundary condition necessary for the 2D-CEM should be available from hðzÞ computed from the all-atom simulation, shown in Figs. 3 and 4. However, due to the ‘‘shape’’ of the inclusion, in this case, a squat cylinder, lipids go over the top of the channel such that hðzÞ is no longer the height of the leaflet in the proper sense (i.e., there is no lipid material below the hydrophobic atom). An alternative method for computing the hydrophobic surface is to compute lipid traces. These traces indicate the shape of the channel by recording the path of the lipid from tail to headgroup. The trace indicates both the proper boundary condition for use by the 3D-CEM and an additional way to compute the slope (from the vector normal to the trace). These traces indicate that h0 ðr0 Þ is positive without the ambiguity of examining hðrÞ. The effective matching condition is confirmed by computing the curvature stresses of the leaflet around gA, quan0 tified by F ð0Þ. This thermodynamic quantity indicates that indeed there is net internal positive curvature of the lipids. 0 The sign of F ð0Þ shows how the free energy would change if the leaflet’s curvature could be individually changed. In 0 this case, a positive value of F ð0Þ means the free energy would go down on a negatively curved leaflet, partially undoing the strain of the internal positive curvature created by the inclusion. Unlike what is seen in the all-atom simulation, the 2D-CEM(s ¼ 0) predicts that there should be minimal curvature stress (see Fig. 5). The 2D-CEM(s free) and 3DCEM values are consistent with the sign of the all-atom simulation, yet the 2D-CEM values are too high, suggesting that the first-shell curvature in the model is too high and that reliance on the apparent agreement of hðrÞ discussed above is unwarranted. In contrast, the 3D-CEM is in quantitative agreement with the all-atom simulation. The final metric for comparison is from the model energetics governing lipid rearrangement around the peptide. In this case, again, the 2D-CEM ðs ¼ 0Þ fails, predicting more enrichment of short (matching) lipids around the gA channel due to the heavy constraint of the slope. The softer energetics of the 2D-CEM(s free) and the 3D-CEM are nicely in line with the observations from the simulation, detailed in Tables 6 and 7. Together, these observations point to the use of the additional parameters of the 3D-CEM to consistently model the interplay between protein shape and membrane deformation energetics. The use of the all-atom simulations to compute curvature stress indicates the matching condition (i.e., a positive slope) without interpreting a molecular feature like hðzÞ at the atomic level. Additionally, by computing local

Biophysical Journal 112, 1198–1213, March 28, 2017 1209

Sodt et al.

lipid directors (lipid ‘‘traces’’) the effective shape of the channel is revealed in a way that can be transferred directly to the 3D-CEM. By combining the 3D-CEM and these simulation observations, the energetics of the membrane surface are determined and validated. The degree of lipid rearrangement expected for a lipid mixture, and thus the predicted lipid enrichment adjacent to the channel, is remarkably similar to that deduced from the all-atom MD simulations in Article I. Yet there is inherent danger in making quantitative conclusions about molecular phenomena using continuum models. In the end, the value of H (or Hd) deduced from experiments is an empirical parameter that, with appropriate care, can be interpreted using a CEM, but there are many uncertainties associated with applying continuum models at the molecular level; previously, all these uncertainties were lumped into the deduced value of s. Questions also remain regarding how far conclusions drawn from the 3D-CEM can be transferred beyond the scope of simulation. In this work, curvature stress of the CEM was validated against molecular dynamics simulations 0 by comparing F ð0Þ computed from lateral stresses. However, in both simulation and CEM, it is possible to extract microscopic information about lateral stress that is unobservable (yet implies observable properties). The finely resolved material parameters of the 3D-CEM, necessary to reproduce the position of the neutral surface and ratio of compression to bending moduli, have a similar ambiguity. A more far-reaching test of the 3D-CEM will be to model interactions between proteins crowded together in heterogeneous bilayers. The bilayer deformation near proteins will be sensitive to the protein-lipid boundary. Thus, coupling between proteins will be similarly affected. The 3D-CEM and lipid-trace technique offer a way to resolve the protein-lipid boundary ambiguity. APPENDIX A: FIRST-SHELL LIPIDS: CURVATURE STRESS Given a point on an axially symmetric surface, fx; y; zg ¼ fr cosðqÞ; r sinðqÞ; hðrÞg, the principal radial and tangential curvatures (cr and ct , respectively) of a single lipid leaflet are (with small-deformation approximations and directions of curvature after the colon)

cr ¼


 z  h00 ðrÞ : 5 1 þ h0 ðrÞÞ2 (33)  3=2 2

h00 ðrÞ 1 þ h0 ðrÞ

and

h0 ðrÞ ct ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiz  r 1 h0 ðrÞ : 5 f sinðqÞ; cosðqÞ; 0g; r 1 þ h0 ðrÞ2 (34) where h0 ðrÞ and h00 ðrÞ are the first and second derivatives of the height. Given hðrÞ, these expressions indicate how to evaluate the curvature modeled by the Helfrich Hamiltonian, and they will be used to relate the

1210 Biophysical Journal 112, 1198–1213, March 28, 2017

curvature stress to the slope of the leaflet where it meets the inclusion, h0 ðr0 Þ. For the upper monolayer, the sign convention has been chosen so that the lipid headgroups are pointing in the positive z direction. First, consider the case of a nearly flat surface with small h0 ðrÞ. In this case, the denominator of cr is ~1 and can be ignored. The direction of this principal curvature is radial and corresponds simply to h00 ðrÞ; thus, this curvature is labeled with ‘‘r.’’ Second, consider the limit where h0 ðrÞ/N, corresponding to a cylinder with lipid headgroups inside. In this case, curvature ct reduces to ð1=rÞ (the curvature of a cylinder), whereas for small h0 ðrÞ, the curvature is proportional to h0 ðrÞ. The direction of ct is perpendicular to the radial curvature, so this curvature is labeled with ‘‘t’’ for tangential. Following (52), in the Monge gauge (for small deviations of hðrÞ from planarity), the total curvature of the system depends on the slope of the boundary condition for where the inclusion meets the leaflet. The total curvature is

V2 hðrÞ ¼ cr þ ct ¼ h00 ðrÞ  h0 ðrÞ=r:

(35)

For a gA channel that is shorter than the bilayer, h00 ðr0 Þ is positive, whereas h0 ðr0 Þ=r0 is negative, defining the sign convention that positive curvature is convex with respect to the headgroups (the reverse of the mathematical sign). Integrating from the channel radius to a distance sufficiently far from the channel (such that deformation from the bulk is negligible) yields, for cr , the radial principal curvature:

Z

2p

N

r0

Z

drrcr ¼ 2p

N

drr½h00 ðrÞ

r0

N

¼ 2p½rh0 ðrÞr0 þ 2p

Z

N

drh0 ðrÞ

r0

¼ 2pr0 h0 ðr0 Þ þ 2p½hðNÞ  hðr0 Þ:

(36)

The result for ct , the tangential principal curvature, is

Z

2p

N

r0

Z

drrct ¼ 2p

N

r0

0 h ðrÞ drr r

¼ 2p½hðNÞ  hðr0 Þ:

(37)

Here, hðNÞ is the height of the bulk, unperturbed lipid. Thus, cr and ct each have components that are proportional to bilayer mismatch but are opposite in sign and equal in magnitude. The radial component has an additional component that is proportional to the slope of the boundary condition. The sum of the two curvature contributions in the Monge gauge is thus

Z

N r0

dr2pr½cr þ ct  ¼ 2pr0 h0 ðr0 Þ

(38)

In the 2D-CEM, the total curvature depends only on the slope of the inclusion-lipid boundary condition, in the limit of small deformations where the Monge gauge applies. Thick bilayers containing lipids with negative spontaneous curvature should have reduced channel lifetime for h0 ðr0 Þ > 0, the case of net positive curvature. That is, lipids that energetically prefer to bend toward their headgroups are forced to bend toward the tails to meet the hydrophobic thickness of the channel. This sensitivity to curvature stress has been observed by Lundbaek et al. (33), in the form of shorter channel lifetimes for systems containing phosphatidylethanolamine lipids. From simulation, this same 0 positive curvature stress is apparent from the change in F ð0Þ (16,17,74), the derivative of the free energy with respect to curvature, upon addition of a mismatched inclusion, as is reported herein.

Modeling the Gramicidin-Lipid Boundary

First-shell lipids: compression stress

6. Brown, M. F. 2012. Curvature forces in membrane lipid-protein interactions. Biochemistry. 51:9782–9795.

The compression energy of lipids in the first shell is modeled as

7. Brown, F. L. H. 2008. Elastic modeling of biomembranes and lipid bilayers. Annu. Rev. Phys. Chem. 59:685–712.

Z

r1



KA;m hðrÞ  h0 2 h0

2 (39)

8. Pan, J., S. Tristram-Nagle, ., J. F. Nagle. 2008. Temperature dependence of structure, bending rigidity, and bilayer interactions of dioleoylphosphatidylcholine bilayers. Biophys. J. 94:117–124.

where r1 is chosen to integrate over the first shell of lipids. Truncating hðrÞ at the first order gives

9. Hamm, M., and M. M. Kozlov. 1998. Tilt model of inverted amphiphilic mesophases. Eur. Phys. J. B. 6:519–528.

(40)

10. Wang, X., and M. Deserno. 2016. Determining the lipid tilt modulus by simulating membrane buckles. J. Phys. Chem. B. 120:6061–6073.

This yields an approximation for Eshell to the first order in h0 ðr0 Þ and Dr:

11. Watson, M. C., E. G. Brandt, ., F. L. H. Brown. 2012. Determining biomembrane bending rigidities from simulations of modest size. Phys. Rev. Lett. 109:028102.

Eshell ¼

r0

2pr

dr;

2 hðrÞ ¼ hp þ h0 ðr0 Þðr  r0 Þ þ O ðr  r0 Þ

Eshell h0 ðr0 ÞDr z1  ¼ 1  ffs ; h0  hp E0

(41)

where r ¼ ðr0 þ r1 Þ=2, Dr ¼ ðr1  r0 Þ, Dh ¼ h0  hp , and



E0 ¼

KA;m Dh h0 2

2 (42)

is the energy density in the immediate vicinity of the inclusion, where compression is the highest. The quantity ffs can be understood as the fraction of the height deformation that is removed within the vicinity of the inclusion. A substantial portion of the compression energy can be removed by an initial membrane deflection from the peptide (this of course requires a curvature penalty).

12. Cantor, R. S. 1999. The influence of membrane lateral pressures on simple geometric models of protein conformational equilibria. Chem. Phys. Lipids. 101:45–56. 13. Botelho, A. V., N. J. Gibson, ., M. F. Brown. 2002. Conformational energetics of rhodopsin modulated by nonlamellar-forming lipids. Biochemistry. 41:6354–6368. 14. Brown, M. F. 1994. Modulation of rhodopsin function by properties of the membrane bilayer. Chem. Phys. Lipids. 73:159–180. 15. Brown, M. F. 1997. Influence of nonlamellar-forming lipids on rhodopsin. In Current Topics in Membranes. R. M. Epand, editor. Academic Press, New York, pp. 285–356. 16. Szleifer, I., D. Kramer, ., S. A. Safran. 1990. Molecular theory of curvature elasticity in surfactant films. J. Chem. Phys. 92:6800–6817. 17. Goetz, R., and R. Lipowsky. 1998. Computer simulations of bilayer membranes: self-assembly and interfacial tension. J. Chem. Phys. 108:7397–7409.

AUTHOR CONTRIBUTIONS

18. Landau, L. D., and E. M. Lifshitz. 1970. Theory of Elasticity. Pergamon Press, Oxford, New York.

A.J.S. and A.H.B. performed the research and wrote the article. All authors (A.J.S., A.H.B., O.S.A., W.I., and R.W.P.) helped design the research and contributed to writing the article.

19. Rand, R. P., N. L. Fuller, ., V. A. Parsegian. 1990. Membrane curvature, lipid segregation, and structural transitions for phospholipids under dual-solvent stress. Biochemistry. 29:76–87.

ACKNOWLEDGMENTS

20. Tian, A., B. R. Capraro, ., T. Baumgart. 2009. Bending stiffness depends on curvature of ternary lipid mixture tubular membranes. Biophys. J. 97:1636–1646.

Computational work was performed on the LoBoS computer cluster of the NHLBI, as well as the Biowulf system of the National Institutes of Health. This work was supported in part by the intramural research program of the National Heart Lung Blood Institute (NHLBI) and the Eunice Kennedy Shriver National Institute of Child Health and Human Development of the National Institutes of Health (NIH), the National Science Foundation (MCB-1157677 and MCB-1727508 to W.I.), and the extramural research program of the National Institutes of Health (GM021342 to O.S.A.).

REFERENCES 1. Sackmann, E. 1984. Physical basis for trigger processes and membrane structures. In Biological Membranes, Vol. 5. D. Chapman, editor.. Academic Press, New York, pp. 105–143. 2. Mouritsen, O. G., and M. Bloom. 1984. Mattress model of lipid-protein interactions in membranes. Biophys. J. 46:141–153. 3. Evans, E., and D. Needham. 1987. Physical properties of surfactant bilayer membranes: thermal transitions, elasticity, rigidity, cohesion and colloidal interactions. J. Phys. Chem. 91:4219–4228. 4. Canham, P. B. 1970. The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell. J. Theor. Biol. 26:61–81. 5. Helfrich, W. 1973. Elastic properties of lipid bilayers: theory and possible experiments. Z. Naturforsch. C. 28:693–703.

21. Patra, M. 2005. Lateral pressure profiles in cholesterol-DPPC bilayers. Eur. Biophys. J. 35:79–88. 22. Siegel, D. P., and M. M. Kozlov. 2004. The gaussian curvature elastic modulus of N-monomethylated dioleoylphosphatidylethanolamine: relevance to membrane fusion and lipid phase behavior. Biophys. J. 87:366–374. 23. Soubias, O., W. E. Teague, Jr., ., K. Gawrisch. 2010. Contribution of membrane elastic energy to rhodopsin function. Biophys. J. 99: 817–824. 24. Sukharev, S., M. Betanzos, ., H. R. Guy. 2001. The gating mechanism of the large mechanosensitive channel MscL. Nature. 409: 720–724. 25. Gullingsrud, J., and K. Schulten. 2004. Lipid bilayer pressure profiles and mechanosensitive channel gating. Biophys. J. 86:3496–3509. 26. Yeagle, P. L., M. Bennett, ., A. Watts. 2007. Transmembrane helices of membrane proteins may flex to satisfy hydrophobic mismatch. Biochim. Biophys. Acta. 1768:530–537. 27. Yin, F., and J. T. Kindt. 2012. Hydrophobic mismatch and lipid sorting near OmpA in mixed bilayers: atomistic and coarse-grained simulations. Biophys. J. 102:2279–2287. 28. Kaiser, H.-J., A. Orlowski, ., K. Simons. 2011. Lateral sorting in model membranes by cholesterol-mediated hydrophobic matching. Proc. Natl. Acad. Sci. USA. 108:16628–16633. 29. Helfrich, P., and E. Jakobsson. 1990. Calculation of deformation energies and conformations in lipid membranes containing gramicidin channels. Biophys. J. 57:1075–1084.

Biophysical Journal 112, 1198–1213, March 28, 2017 1211

Sodt et al. 30. Huang, H. W. 1986. Deformation free energy of bilayer membrane and its effect on gramicidin channel lifetime. Biophys. J. 50:1061–1070. 31. Nielsen, C., M. Goulian, and O. S. Andersen. 1998. Energetics of inclusion-induced bilayer deformations. Biophys. J. 74:1966–1983. 32. Sodt, A. J., and R. W. Pastor. 2014. Molecular modeling of lipid membrane curvature induction by a peptide: more than simply shape. Biophys. J. 106:1958–1969. 33. Lundbaek, J. A., A. M. Maer, and O. S. Andersen. 1997. Lipid bilayer electrostatic energy, curvature stress, and assembly of gramicidin channels. Biochemistry. 36:5695–5701. 34. Lundbaek, J. A., S. A. Collingwood, ., O. S. Andersen. 2010. Lipid bilayer regulation of membrane protein function: gramicidin channels as molecular force probes. J. R. Soc. Interface. 7:373–395. 35. Miloshevsky, G. V., and P. C. Jordan. 2004. Gating gramicidin channels in lipid bilayers: reaction coordinates and the mechanism of dissociation. Biophys. J. 86:92–104.

52. Ring, A. 1996. Gramicidin channel-induced lipid membrane deformation energy: influence of chain length and boundary conditions. Biophys. Biochim. Acta. 1278:147–159. 53. Kim, T., K. I. Lee, ., W. Im. 2012. Influence of hydrophobic mismatch on structures and dynamics of gramicidin a and lipid bilayers. Biophys. J. 102:1551–1560. 54. Lee, K. I., R. W. Pastor, ., W. Im. 2013. Assessing smectic liquidcrystal continuum models for elastic bilayer deformations. Chem. Phys. Lipids. 169:19–26. 55. Nelder, J. A., and R. Mead. 1965. A simplex method for function minimization. Comput. J. 7:308–313. 56. Harroun, T. A., W. T. Heller, ., H. W. Huang. 1999. Experimental evidence for hydrophobic matching and membrane-mediated interactions in lipid bilayers containing gramicidin. Biophys. J. 76:937–945. 57. Campelo, F., H. T. McMahon, and M. M. Kozlov. 2008. The hydrophobic insertion mechanism of membrane curvature generation by proteins. Biophys. J. 95:2325–2339.

36. Mondal, S., G. Khelashvili, ., H. Weinstein. 2011. Quantitative modeling of membrane deformations by multihelical membrane proteins: application to G-protein coupled receptors. Biophys. J. 101: 2092–2101.

58. Kozlov, M. M. and M. Winterhalter. Elastic moduli for strongly curved monoplayers. Position of the neutral surface. J. Phys. II. 1: 1077–1084 .

37. Elliot, J. R., D. Needham, ., D. A. Haydon. 1983. The effects of bilayer thickness and tension on gramicidin single-channel lifetime. Biophys. Biochim. Acta. 735:95–103.

59. Prak, D. J. L., E. K. Brown, and P. C. Trulove. 2013. Density, viscosity, speed of sound, and bulk modulus of methyl alkanes, dimethyl alkanes, and hydrotreated renewable fuels. J. Chem. Eng. Data. 58:2065–2075.

38. Lundbaek, J. A., and O. S. Andersen. 1999. Spring constants for channel-induced lipid bilayer deformations. Estimates using gramicidin channels. Biophys. J. 76:889–895.

60. May, S., Y. Kozlovsky, and M. M. Kozlov. 2004. Tilt modulus of a lipid monolayer. Eur. Phys. J. E. 14:299–308.

39. Klauda, J. B., R. M. Venable, ., R. W. Pastor. 2010. Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types. J. Phys. Chem. B. 114:7830–7843. 40. Venable, R. M., F. L. H. Brown, and R. W. Pastor. 2015. Mechanical properties of lipid bilayers from molecular dynamics simulation. Chem. Phys. Lipids. 192:60–74. 41. Botan, A., F. Favela-Rosales, ., J. Tynkkynen. 2015. Toward atomistic resolution structure of phosphatidylcholine headgroup and glycerol backbone at different ambient conditions. J. Phys. Chem. B. 119: 15075–15088. 42. Yoo, J., and Q. Cui. 2013. Three-dimensional stress field around a membrane protein: atomistic and coarse-grained simulation analysis of gramicidin A. Biophys. J. 104:117–127. 43. Argudo, D., N. P. Bethel, ., M. Grabe. 2016. Continuum descriptions of membranes and their interaction with proteins: towards chemically accurate models. Biochim. Biophys. Acta. 1858 (7 Pt. B):1619–1634. 44. Nielsen, C., and O. S. Andersen. 2000. Inclusion-induced bilayer deformations: effects of monolayer equilibrium curvature. Biophys. J. 79:2583–2604. 45. Brannigan, G., and F. L. H. Brown. 2007. Contributions of Gaussian curvature and nonconstant lipid volume to protein deformation of lipid bilayers. Biophys. J. 92:864–876. 46. Safran, S. A. 1994. Statistical Thermodynamics of Surfaces, Interfaces, and Membranes. Westview Press, Boulder, CO. 47. Rawicz, W., K. C. Olbrich, ., E. Evans. 2000. Effect of chain length and unsaturation on elasticity of lipid bilayers. Biophys. J. 79:328–339. 48. Evans, E., W. Rawicz, and B. A. Smith. 2013. Back to the future: mechanics and thermodynamics of lipid biomembranes. Faraday Discuss. 161:591–611. 49. Henriksen, J. R., and J. H. Ipsen. 2004. Measurement of membrane elasticity by micro-pipette aspiration. Eur Phys J E Soft Matter. 14:149–167. 50. Dimova, R. 2014. Recent developments in the field of bending rigidity measurements on membranes. Adv. Colloid Interface Sci. 208: 225–234. 51. Bouvrais, H., L. Duelund, and J. H. Ipsen. 2014. Buffers affect the bending rigidity of model lipid membranes. Langmuir. 30:13–16.

1212 Biophysical Journal 112, 1198–1213, March 28, 2017

61. May, S. 2000. Protein-induced bilayer deformations: the lipid tilt degree of freedom. Eur. Biophys. J. 29:17–28. 62. Nocedal, J. 1980. Updating quasi-Newton matrices with limited storage. Math. Comput. 35:773–782. 63. Jo, S., T. Kim, ., W. Im. 2008. CHARMM-GUI: a web-based graphical user interface for CHARMM. J. Comput. Chem. 29:1859–1865. 64. Lee, J., X. Cheng, ., W. Im. 2016. CHARMM-GUI input generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM simulations using the CHARMM36 additive force field. J. Chem. Theory Comput. 12:405–413. 65. Wu, E. L., X. Cheng, ., W. Im. 2014. CHARMM-GUI Membrane Builder toward realistic biological membrane simulations. J. Comp. Chem. 35:1997–2004. 66. Townsley, L. E., W. A. Tucker, ., J. F. Hinton. 2001. Structures of gramicidins A, B, and C incorporated into sodium dodecyl sulfate micelles. Biochemistry. 40:11676–11686. 67. Mackerell, A. D., Jr., M. Feig, and C. L. Brooks, 3rd. 2004. Extending the treatment of backbone energetics in protein force fields: limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J. Comput. Chem. 25:1400–1415. 68. Andersen, H. C. 1983. RATTLE: A ‘‘velocity’’ version of the SHAKE algorithm for molecular dynamics calculations. J. Comput. Phys. 52:24–34. 69. Darden, T., D. York, and L. Pedersen. 1993. Particle mesh Ewald: an N$log(N) method for Ewald sums in large systems. J. Chem. Phys. 98:10089–10092. 70. Brooks, B. R., C. L. Brooks, 3rd, ., M. Karplus. 2009. CHARMM: the biomolecular simulation program. J. Comput. Chem. 30:1545–1614. 71. Lindahl, E., and O. Edholm. 2000. Spatial and energetic-entropic decomposition of surface tension in lipid bilayers from molecular dynamics simulations. J. Chem. Phys. 113:3882–3893. 72. Sonne, J., F. Y. Hansen, and G. H. Peters. 2005. Methodological problems in pressure profile calculations for lipid bilayers. J. Chem. Phys. 122:124903. 73. Ollila, O. H. S., H. J. Risselada, ., S. J. Marrink. 2009. 3D pressure field in lipid membranes and membrane-protein complexes. Phys. Rev. Lett. 102:078101.

Modeling the Gramicidin-Lipid Boundary 74. Sodt, A. J., and R. W. Pastor. 2013. Bending free energy from simulation: correspondence of planar and inverse hexagonal lipid phases. Biophys. J. 104:2202–2211. 75. Metropolis, N., A. W. Rosenbluth, ., E. Teller. 1953. Equation of state calculations by fast computing machines. J. Chem. Phys. 21:1087– 1092. 76. Partenskii, M. B., and P. C. Jordan. 2002. Membrane deformation and the elastic energy of insertion: perturbation of membrane elastic constants due to peptide insertion. J. Chem. Phys. 117:10768–10776.

77. Greisen, P., Jr., K. Lum, ., J. Lundbæk. 2011. Linear rate-equilibrium relations arising from ion channel-bilayer energetic coupling. Proc. Natl. Acad. Sci. USA. 108:12717–12722. 78. Nagle, J. F., M. S. Jablin, ., K. Akabori. 2015. What are the true values of the bending modulus of simple lipid bilayers? Chem. Phys. Lipids. 185:3–10. 79. Siegel, D. P. 2008. The Gaussian curvature elastic energy of intermediates in membrane fusion. Biophys. J. 95:5200–5215.

Biophysical Journal 112, 1198–1213, March 28, 2017 1213