This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.105.075838v1 90/5/1501    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Brannigan, G. Articles by Brown, F. L. H. Search for Related Content PubMed PubMed Citation Articles by Brannigan, G. Articles by Brown, F. L. H.

A Consistent Model for Thermal Fluctuations and Protein-Induced Deformations in Lipid Bilayers

Grace Brannigan ^* and Frank L. H. Brown 

^* Department of Physics and Astronomy, and Department of Chemistry and Biochemistry, University of California, Santa Barbara, California

Correspondence: Address reprint requests to F. L. H. Brown, Tel.: 805-893-5494; Fax: 805-893-4120; E-mail: fbrown{at}chem.ucsb.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION A GENERAL MODEL FLUCTUATION SPECTRA OF... COMPARISON TO SIMULATION DATA PROTEIN-INDUCED DEFORMATION... MOLECULAR SIMULATION DISCUSSION APPENDIX A: APPROXIMATIONS TO... APPENDIX B: CURVE-FITTING AND... ACKNOWLEDGEMENTS REFERENCES

 
We present an elastic Hamiltonian for membrane energetics that captures bilayer undulation and peristaltic deformations over all wavelengths, including the short wavelength protrusion regime. The model implies continuous functional forms for thermal undulation and peristaltic amplitudes as a function of wavelength and predicts previously overlooked relationships between these curves. Undulation and peristaltic spectra display excellent agreement with data from both atomistic and coarse-grained models over all simulated length scales. Additionally, the model accurately predicts the bilayer's response to a cylindrical protein inclusion as observed in coarse-grained simulation. This elastic response provides an explanation for gramicidin ion channel lifetime versus membrane thickness data that requires no fit constants. The physical parameters inherent to this picture may be expressed in terms of familiar material properties associated with lipid monolayers. Inclusion of a finite monolayer spontaneous curvature is essential to obtain fully consistent agreement between theory and the full range of available simulation/experimental data.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION A GENERAL MODEL FLUCTUATION SPECTRA OF... COMPARISON TO SIMULATION DATA PROTEIN-INDUCED DEFORMATION... MOLECULAR SIMULATION DISCUSSION APPENDIX A: APPROXIMATIONS TO... APPENDIX B: CURVE-FITTING AND... ACKNOWLEDGEMENTS REFERENCES

 
At and near physiological temperatures, lipid bilayers exhibit significant thermal fluctuations in microscopic structure, composition, and shape as dictated by equilibrium statistical mechanics (1,2). Membranes are not static, flat homogeneous structures—not only because of metabolic activity and biological structures at the plasma membrane surface (cytoskeleton, caveolae, lipid rafts, coated pits, etc.), but also because of these purely physical considerations. Although living cells are certainly not equilibrium structures, it is important to fully understand the thermal behavior of model membrane systems as a preliminary step toward unraveling biologically relevant phenomena at membrane surfaces.

Thermal fluctuations in lipid bilayers have been implicated in a variety of biophysical phenomena, including (but not limited to) steric repulsions between proximal bilayers (3,4), shape fluctuations of the red blood cell (5), cellular motility (6), and entropically driven interactions between integral membrane proteins (7,8). Traditionally, our theoretical understanding of such phenomena has rested upon simplified analytical theories, in the spirit of work by Helfrich (9) and others (3,10,11). Since these theories describe the bilayer by one or more continuous fields in space without any atomic/molecular level resolution, we will refer to them as elastic pictures. Physical properties needed in the formulation of such theories (elastic moduli, interfacial tensions, etc.) are typically guessed or inferred (often indirectly) from experiment. More recently, molecular dynamics simulation has evolved as a potential means to connect theoretical models with specific lipid bilayers by providing physical parameters directly from computer experiments.

In principle, it should be possible to extract the parameters inherent to elastic models directly from atomistic simulations. In such a scheme, differences in behavior due to lipid composition would be fully predicted by detailed simulations of chemically distinct bilayers over relatively small length scales. In practice, the correspondence between elastic models and molecular simulations has yet to be fully developed. Although the extraction of bilayer bending moduli and surface tensions through analysis of thermal membrane height fluctuations (undulations) is well established (12–19), the corresponding interpretation of bilayer thickness (peristaltic) fluctuations is less common (13–15). (See Fig. 1 for an explanation of height versus thickness fluctuations.) Perhaps one reason for this is that current protocols for fitting peristaltic data involve separately analyzing simulations over multiple wavelength regimes, without clear rules for how to effect such a separation. Furthermore, one expects that the elastic properties involved in both peristaltic and height fluctuations derive from the same source, yet present theories do not fully account for the interrelations between these two types of deformation. Finally, it should be noted that there is some discrepancy in how various groups treat data stemming from the short wavelength protrusion regime (Fig. 1 A). Although numerous theoretical studies have considered long-wavelength undulation modes, long-wavelength peristaltic modes, and protrusion modes in separate contexts (3,9–11,13), to date no single unifying theory has been advanced that can account for all of these behaviors. So, although the physical basis for these behaviors is fairly well established, it is unclear how to best interpret data quantitatively, since it is not fully clear how all of these fluctuation modes are coupled and it is not known how thermal undulation and peristaltic amplitudes are expected to relate to one another.

View larger version (18K): [in this window] [in a new window]   FIGURE 1  (A) Fluctuation modes of lipid bilayers (3,11). Over long lengthscales, fluctuations are dominated by mesoscopic deformations involving the concerted motion of many lipids. At wavelengths of several nanometers and longer, molecular-level details are unimportant, with fluctuations dominated by these long-wavelength or bending modes. Molecular-level roughness is, of course, unavoidable at sufficiently short wavelengths. The bilayer motions associated with molecular level fluctuations are traditionally denoted "protrusions". The two leaflets need not move in phase, which can result in a nonuniform thickness of the bilayer in both the protrusion and bending regimes. (B) Definition of height and thickness fluctuations. Shape fluctuations of the bilayer are conveniently decomposed into height and thickness contributions. In this work, we adopt the convention that the membrane is fluctuating around an average flat configuration with normal in the z direction. We denote the midplane between monolayer leaflets as a function of x,y coordinates as the height field of the membrane (denoted h, see Eq. 5). The distance between monolayer leaflets as a function of x,y position is referred to as the "thickness field of the membrane" (denoted 2t + 2t0, the fluctuating variable t is defined as the thickness of a monolayer measured relative to the average tensionless thickness to conform with convention; see Eq. 6). In practice, the position of the individual monolayers is dictated by specific atomic groups associated with the interface between polar and hydrophobic groups along the lipid chain. We shall use the terms "height fluctuation" and "undulation" interchangeably to describe deviations of the height field away from the flat reference state. Likewise, we shall refer to "thickness fluctuations" and "peristaltic modes" interchangeably. In contrast to some authors, we allow these terms to refer to fluctuations over all wavelengths—i.e., bending and protrusions both contribute to undulations and peristaltic modes.

View larger version (30K): [in this window] [in a new window]   FIGURE 5  Inclusion-induced deformation. A symmetric transmembrane protein with hydrophobic residues around its periphery and a thickness exceeding that of the surrounding membrane will tend to distort the bilayer as shown (not to scale). A protein thinner than the surrounding membrane is expected to induce the opposite effect. The nonmonotonic healing of the membrane thickness back to equilibrium separation as one moves away from the inclusion is predicted theoretically by our model (which nearly reduces to the theory of Aranda-Espinoza et al. (31) in the context of inclusion-induced deformations), and observed in our simulations and similar coarse-grained simulations (36,37).

Other applications are certainly possible as well. Although many of the underlying components of our theory have been discussed in specific contexts previously, the unified formulation we present is new and appealing in its ability to consistently subsume a variety of different physical phenomena. This picture naturally resolves several open questions and inconsistencies, while serving as a convenient means to analyze simulated and experimental data. The picture we advance is fully consistent with an array of such data. Specifically, we call attention to the following aspects of this work:

1. We derive an expression for thermal peristaltic fluctuations from an underlying elastic model. The corresponding expressions for thermal undulations are derived from the same model and correspond to the expected (11) Helfrich behavior at long wavelengths and protrusion behavior over molecular wavelengths. Since both thickness and height fluctuations are derived from the same starting point, we find explicit correspondence between physical parameters quantifying both types of fluctuation. In particular, we predict identical bending moduli and protrusion-associated constants for both phenomena.
   
2. The possibility of nonvanishing monolayer spontaneous curvature is central to our model. Under certain parameter regimes, this leads to nonmonotonic behavior for thermal peristaltic fluctuation amplitudes as a function of wavelength (undulation modes are not affected). Such behavior has been observed in fully atomic simulations of dipalmitoylphosphatidylcholine (DPPC) (13) and sphingomyelin (SM) (15) bilayers, but has previously been attributed to poor sampling. Our treatment provides an appealing alternative explanation.
   
3. We have used our theory to fit existing fully atomic MD simulation fluctuation data from the literature for three different lipid bilayers (DPPC (13), glycerolmonoolein (GMO) (14), and SM (15)) and a coarse-grained (CG) bilayer model developed in our group (16). There are a total of five fit parameters involved in this process, with clearly defined physical meanings. Previous fitting schemes have all involved at least this many constants, yet did not account for the important contribution of finite monolayer spontaneous curvature. The obtained numerical values for these constants are physically reasonable and the quality of our fits to the data sets are universally very good (see Fig. 4).
   
4. We critically analyze the results of a coarse-grained simulation of a cylindrical protein embedded within an otherwise homogeneous bilayer. The proposed elastic picture quantitatively predicts the bilayer's response to protein-induced deformation, using elastic constants obtained (within error bars) via analysis of thermal fluctuations in the homogeneous bilayer. In other words, the elastic picture we advance finds fully consistent agreement between two very different simulations.
   
5. Using elastic properties derived from thermal fluctuation data of GMO, we predict the effect of bilayer thickness on gramicidin-A ion-channel lifetimes for monoglyceride-based bilayers. Our procedure involves no fit parameters and our model predicts relative lifetimes in good agreement with experiment.
   
6. Applied to protein-induced deformations, our model nearly reduces to the picture advanced by Aranda-Espinoza et al. (31). In contrast to that picture, we include the possibility of microscopic protrusions and solve for the thermal average of the deformation profile (numerically) as opposed to identifying the elastic minimum (analytically). These differences are expected to produce negligible effects for three of the four bilayers analyzed in this study and are explicitly shown to produce negligible effect for our CG simulation model. In this sense, we have provided the first direct validation of Aranda-Espinoza's theory for predicting membrane shape around symmetric inclusions with hydrophobic mismatch. In the case of GMO, protrusion-bending coupling is relatively strong and we predict that protrusions will affect the deformation profile, leading to quantitative disagreement with Aranda-Espinoza et al. (31). In the absence of appropriate simulation data on GMO, this prediction remains unverified.
   
7. Applied to thermal membrane fluctuations, our model predicts continuous functional forms for both undulation and peristaltic amplitudes over all wavelengths. This resolves the practical shortcomings of previously introduced piecewise fitting techniques and, as noted above, explains certain interrelations between these data sets.
   

View larger version (24K): [in this window] [in a new window]   FIGURE 4  Height (|hq|2) and thickness (|tq|2) fluctuations for DPPC (13), GMO (14), SM (15), and a coarse-grained model (CG). Simulation data are displayed as circles. The lines represent best fits of the data to the expressions in Eq. 32. Fit parameters are in Table 1.

	A GENERAL MODEL

TOP ABSTRACT INTRODUCTION A GENERAL MODEL FLUCTUATION SPECTRA OF... COMPARISON TO SIMULATION DATA PROTEIN-INDUCED DEFORMATION... MOLECULAR SIMULATION DISCUSSION APPENDIX A: APPROXIMATIONS TO... APPENDIX B: CURVE-FITTING AND... ACKNOWLEDGEMENTS REFERENCES

 
We derive here our model for bilayer deformations. Because we seek to explain both height and thickness deformations (Fig. 1), our considerations begin with a bilayer composed of two opposing coupled monolayers. The two leaflets do not necessarily bend in unison, which leads to both height and thickness fluctuations/deformations. Our starting point is the general theory for surfactant monolayers presented in Safran (1), but we retain fluctuations in area per lipid. We assume a constant volume condition for hydrocarbon tails and additionally assume that lipids across from one another in opposing leaflets share the same local area/lipid. This leads to a theory for membrane elasticity in which thickness and height deformations are uncoupled. The thickness deformations obey energetics consistent with the picture developed in Aranda-Espinoza et al. (31), whereas height fluctuations are consistent with standard Helfrich energetics. We extend this approach to include microscopic protrusions as in the theory of Lipowsky and Grotehans (11), but the formulation discussed here is more general since it includes the contribution of peristaltic bending modes.

Each monolayer can exhibit bending deformations described by fields z^{(1, 2)} (denoted as such because we always assume the x,y plane as the reference configuration for bilayer midplane), and is also subject to microscopic noise, or protrusions, described by fields ^{(1, 2)}. By convention, we always take the top leaflet to be monolayer (1). Consequently the bilayer is described by four separate (but coupled) fields: z⁽¹⁾(x,y), z⁽²⁾(x,y), ⁽¹⁾(x,y), and ⁽²⁾(x,y) (Fig. 2). From this point on and for notational simplicity, we assume the x,y dependence without explicitly writing it. It is convenient to define

(1)

(2)

(3)

(4)

Then, as commonly measured in simulations, the height h of the bilayer midplane and deviations in the bilayer thickness 2t are described by

(5)

(6)

respectively. From a theoretical perspective, we imagine these fields to reflect precisely defined interfaces between lipid hydrocarbon tails and surrounding water (Fig. 2). From a practical perspective, these fields are extracted from simulation by triangulating a surface using atoms along the lipid chain to represent the position of such an interface (see Fluctuation Spectra of Homogeneous Membranes for further elaboration). Note that our thickness field t refers to fluctuations of half the bilayer thickness and that these fluctuations are measured relative to t₀, the half-bilayer thickness for a flat sheet in its minimal energy configuration at vanishing tension. This seemingly odd definition has been adopted to facilitate later connection with thickness fluctuation spectra, as previously reported in the literature.

View larger version (23K): [in this window] [in a new window]   FIGURE 2  Defining the elastic model for bilayer deformations. Microscopic fluctuations (1, 2)(x,y) roughen the molecularly smooth interface z(1, 2)(x,y) between each leaflet and water. The values of (1, 2) are defined relative to z(1, 2), respectively. Note that it is the sum z(1, 2) + (1, 2) that defines the interface between polar and hydrophobic groups for each leaflet, as discussed in Fig. 1. The fields z(1, 2) should be regarded as mesoscopic fields reflecting an implicit local averaging over molecular fluctuations in the sense of Landau order parameters. At short wavelengths there are too few molecules to provide a smooth coarse-grained average, and we include the effect of these molecular fluctuations via the fields (1, 2). This is essentially the picture adopted in Lipowsky and Grotehans (11).

f

Bending contribution
To treat bending energetics, we temporarily neglect the fields and treat the bilayer as the two opposing elastic sheets z⁽¹⁾ and z⁽²⁾. As in Safran (1), we Taylor-expand the free energies per molecule to quadratic order in mean curvature H and molecular area deviation ( – ₀),

(7)

The two monolayers share identical material properties; however, the bottom leaflet (2) is inverted relative to the orientation of the top (1), which accounts for the sign differences in terms with first-order curvature contributions. The value is the molecular free energy for a flat monolayer evaluated at the area per molecule ₀ that minimizes the free energy under conditions of vanishing applied tension. The primes represent derivatives with respect to evaluated at ₀. Subscripts simply reflect the power of H that the constants precede. Note that ⁽¹⁾ and ⁽²⁾ are the areas per molecule as measured perpendicular to the local monolayer normal—i.e.,

(8)

(9)

where _xy is area per molecule projected onto the reference x,y plane. As discussed below, we assume this area to be locally identical for the two leaflets.

We bind the two monolayers by requiring conservation of volume for hydrophobic lipid tails and assuming that lipids directly beneath one another (i.e., same x,y coordinates) in the bilayer configuration share the same projected area/molecule (and hence the same thickness by virtue of volume conservation). Stated mathematically, we assume an equation of state for lipids relating local monolayer thickness to local molecular area to be

(10)

Note that we have made no distinction in thickness or molecular area between the two leaflets, since we assume these quantities to be locally identical for the two monolayers. We comment that Eq. 10 represents only approximate conservation of volume since we do not include the effect of protrusions in this expression and we have neglected contributions due to surface slope and curvature. Furthermore, fixing identical local thickness between opposing monolayers represents a seemingly harsh constraint. The scheme we adopt has the advantage of mathematical simplicity, both in formulation and final results. These assumptions naturally lead to a decoupling between peristaltic and undulation modes with identical bending moduli characterizing these two types of deformation. Similar treatments invoking different assumptions at this stage predict different bending moduli and/or coupling between height and thickness fluctuations (unpublished work); such effects do not seem to be supported by simulation or experiment. In any event, the ultimate justification of these simplifications is the correspondence we find between simulation data and the predictions of our theoretical model. The following sections will demonstrate the correspondence to be very good.

Implementation of Eq. 10 on the expressions in Eq. 7 is handled conveniently in the Monge gauge. The bilayer free energy per area projected onto the x,y plane, f_z, is calculated by noting, for the total bending energy of the bilayer,

(11)

where A is the area of the surfaces projected onto the x,y plane. Implementing the Monge representation for small curvatures (()) in our monolayer free energy expressions, and truncating all expansions to second-order in deformation fields, we arrive at

(12)

It is clear from this expression that undulations and peristaltic excitations are decoupled. We identify one of the constants appearing in the above expression, by comparing the undulation term (the one involving z⁺) with the usual Helfrich formula (1,9). The term associated with involves compression of the bilayer, since z^– and ( – ₀) are connected through volume conservation, which allows us to identify this constant with the usual bilayer compression modulus. The term in braces on the first line integrates to the total number of lipids in the bilayer, and so we identify with the lipid chemical potential. Under conditions of vanishing tension, this quantity vanishes. The remaining constants are clearly related to the spontaneous curvature of the component monolayers, as can be seen by examining Eq. 7 (1,31). Explicitly, we deduce

(13)

For future convenience, the values k_A and k_c are defined here as the compressibility modulus and bending modulus for the bilayer. The analogous quantities associated with the monolayers are obtained by dividing the bilayer values in half. The spontaneous curvature c₀ and area derivative of the spontaneous curvature c'₀ c₀/ are the values associated with the individual monolayers. As is evident from Eq. 12, the derived behavior of the bilayer for two identical opposing leaflets yields vanishing bilayer spontaneous curvature regardless of the monolayer values. Spontaneous curvature of the monolayers manifests itself only through peristaltic deformations and fluctuations. Our convention for sign of c₀ insures that positive values encourage a monolayer to form micelles, and that negative values favor reverse micelles.

Adopting the notation specified above, Eq. 12 becomes

(14)

where we have defined

(15)

The pieces of Eq. 14 are easily interpreted. The first term consists of the usual Helfrich bending component associated with a tensionless fluid surface with vanishing spontaneous curvature. It is this portion of the energetics that is usually considered in long wavelength studies of membranes, where peristaltic fluctuations are assumed to be unimportant. The remaining terms describe contributions to the energetics due to out-of-phase motion (peristaltic deformations) of the two monolayers. The second term corresponds to the energetic cost of area stretching/compression that accompanies a thickness fluctuation. The remaining terms reflect bending energetics due to the peristaltic modes analogous to the undulation expression. Note, however, the terms resembling (thickness-dependent) nonvanishing spontaneous curvature contributions for peristaltic deformations. Spontaneous curvature terms vanish in bilayer undulations (formed from identical monolayers), since the contributions from opposing leaflets exactly cancel. The opposite is true for peristaltic fluctuations—the two monolayers' energetics reinforce one another, and we expect to see a contribution to the peristaltic modes reflecting this fact (unless the monolayers themselves have a vanishing spontaneous curvature and vanishing area derivative of c₀).

Protrusion contribution
Thus far, we have neglected the microscopic protrusion fields ^(1,2). We expect the interface between monolayers and water to be subject to microscopic noise, not represented in the coarse-grained bending energetics discussed in the previous section. At the oil-water interface, this noise has the potential to affect membrane energetics by altering the water structuring proximal to the membrane. Although more refined approaches are possible (11), we adopt an elastic description of this noise similar to that proposed by Lipowsky and Grotehans (11). The free energy per unit area is composed of a surface area term, reflecting the interfacial tension between hydrocarbon tails and water, plus a binding term to keep the fluctuations localized to the coarse-grained average fields z^(1,2),

(16)

(17)

where the additional surface area introduced by the fields is, to quadratic order,

(18)

(19)

We subtract off the areas of the bare sheets' z^(1,2) contribution because the interfacial free-energy-associated molecularly smooth shape changes are incorporated within the free energy f_z. In other words, if the membrane is microscopically smooth, this means = 0 everywhere, and f should vanish.

The interfaces are assumed not to affect one another beyond those contributions seen in f_z, so the total contribution of protrusion energetics is the sum from both interfaces: f = f⁽¹⁾ + f⁽²⁾. In terms of our symmetric/antisymmetric variables (Eqs. 1 and 2), the protrusion contribution to the free energy is

(20)

The main difference between Eq. 20 and the expression of Lipowsky and Grotehans is the presence of the z^– · ^– term. Lipowsky and Grotehans treated the bilayer as two microscopic fields grafted to one central elastic sheet, so z^– (the distance between elastic sheets) was necessarily zero everywhere. An additional minor difference is that in our model, ⁺ and ^– are naturally bound by the same constant (k), whereas Lipowsky and Grotehans allow for two different constants. Given the present formulation, the use of a single k constant is physically required. In some sense, the idea of two different protrusion constants comes about in our picture naturally through the incorporation of peristaltic bending modes (z^–). The relation between peristaltic modes and undulations at short wavelengths is predicted in our model, based on the bending elastic constants previously discussed and two constants associated with protrusions.

The total free energy density for the membrane is f = (f_z + f), yielding the entire free energy for the bilayer in a given configuration as

(21)

We emphasize again that A reflects the area of the lipid bilayer projected to the x,y plane. If the bilayer encompasses regions devoid of lipids (as when a protein is embedded), the integral reflects that portion of space occupied by the lipids. Terms in the first line (+ terms) affect height fluctuations, whereas terms in the second and third lines (– terms) affect thickness fluctuations and the two portions are entirely decoupled. This expression is one of the main results of this article. All other calculations follow from this expression for bilayer energetics. We can immediately apply this Hamiltonian to both thermal fluctuation spectra for bilayers (next section) and the deformation profiles induced by membrane proteins (see Protein-Induced Deformation Profiles). For the most part, individual terms appearing within Eq. 21 (or slight variations) have appeared in one or more previous theoretical treatments of bilayer and/or monolayer energetics. We believe this to be the first time all relevant physical effects have been encompassed within a single framework to describe both undulation and peristaltic excitations over all wavelengths, including the protrusion regime. This formulation has significant practical and theoretical advantages, to be elaborated upon in the following sections.

	FLUCTUATION SPECTRA OF HOMOGENEOUS MEMBRANES

TOP ABSTRACT INTRODUCTION A GENERAL MODEL FLUCTUATION SPECTRA OF... COMPARISON TO SIMULATION DATA PROTEIN-INDUCED DEFORMATION... MOLECULAR SIMULATION DISCUSSION APPENDIX A: APPROXIMATIONS TO... APPENDIX B: CURVE-FITTING AND... ACKNOWLEDGEMENTS REFERENCES

 
The fluctuation spectra characterize the thermal height fluctuations and thickness fluctuations of the bilayer as a function of wavevector q. Measurement of the height fluctuations via simulation has become a standard (12–19) method to calculate the effective bending rigidity and surface tension of homogeneous fluid membranes. Fewer studies (13–15) on thickness fluctuations have been conducted, perhaps because there are few direct experimental counterparts. The present model predicts continuous functional forms for both and and represents an immediate application of our Hamiltonian. The following subsections derive the exact predicted forms for the two spectra and go on to suggest a physically motivated simplification that makes the expressions more compact and simpler to interpret. We demonstrate the close agreement between the full and approximate forms for reasonable physical constants in Appendix A. Using the approximate forms (exact forms yield identical results to within error bars), we fit four sets of simulation data from different lipid bilayer simulations. The calculated fits display excellent agreement with simulation results and predict physically reasonable elastic constants.

Theory
To derive the fluctuation spectrum for a homogeneous membrane, we consider the Hamiltonian in Eq. 21 for a membrane in a square periodic box of area A = L². We use the conventions for the Fourier-transform pair of

(22)

(23)

for an arbitrary function g and q = 2(n, m)/L for the integers n, m = 0, ±1, ...N, where N is dictated by a short wavelength cutoff. In Fourier space, the Hamiltonian equation (Eq. 21) becomes

(24)

Note that since

(25)

under periodic boundary conditions, one term from the general expression does not appear in Eq. 24. The remaining expression implies that one cannot measure c₀ and c₀' independently using the fluctuation spectra, but only the linear combination contained within . To measure c₀ independently, we use the membrane stress profile (Fig. 3).

View larger version (13K): [in this window] [in a new window]   FIGURE 3  Stress profile for homogeneous coarse-grained bilayer with 128 lipids at constant vanishing tension. The monolayer spontaneous curvature is related to (z) through Eq. 33.

(26)

It is well established (12–19) that for long wavelengths (small q), the height fluctuations of a nearly flat bilayer at zero tension follow

(27)

where is the effective bending rigidity, essentially defined by Eq. 27. Expanding Eq. 26 around small q,

(28)

we recover the form of Eq. 27, but with a renormalized bending rigidity,

(29)

This relationship demonstrates that, as first derived by Lipowsky and Grotehans (11), protrusions lower the effective bending rigidity relative to the bare value. Note that remains positive only if 2²/(kk_c) < 1, which sets obvious limitations on the relative magnitudes of the elastic constants as discussed by Lipowsky and Grotehans (11). The collection of constants 2²/(kk_c) represents a measure of the coupling strength between protrusions and bending modes. If this coupling were to become very strong, our entire physical picture would be suspect. We expect (and find in Comparison to Simulation Data, below) this number to be quite small for actual lipid systems.

Decoupled protrusion/bending approximation
An obvious way to interpret fluctuation spectra would be to fit data to the full forms of Eq. 26. This is possible and we include such curves in Appendix A. The algebraically complex nature of these expressions makes interpretation difficult, however—and it would be desirable to have simpler approximate expressions. Previous analyses of simulated height fluctuation spectra have involved functions that display the same small q and large q behavior as in Eq. 26, such as (13–15)

(30)

where q_c is typically on the order of the membrane thickness, or (12)

(31)

The limiting forms suggested above are easily predicted from Eq. 26 in the two extreme regimes of q and q 0. Given actual lipid systems, however, these expressions do not faithfully reproduce the results of full calculations. In Appendix A, we compare Eqs. 31 and 26 for several parameter sets, and demonstrate that Eq. 31 is not a dependable approximation for wavelengths on the order of the membrane thickness.

We suggest alternate approximate forms for the spectra based on the assumption of decoupling between protrusion and bending modes. When all terms containing and are set to zero in Eq. 24, the spectra become

(32)

The suitability of these expressions for describing numerous bilayer systems is demonstrated in Comparison to Simulation Data, below, and their numeric equivalence to the expressions in Eq. 26 for a range of constants is demonstrated in Appendix A. We note that these expressions could have been written down immediately by realizing that h and t both represent the sum of protrusion and bending fields (Eqs. 5 and 6). If these fields are independent of one another (i.e., no coupling), the variance of these sums is the sum of the variances associated with protrusion fluctuations and bending fluctuations. The two terms in both of the expressions in Eq. 32 reflect individual contributions from bending and protrusions. The most novel of all these terms is the bending contribution of the thickness spectra. To our knowledge, the incorporation of spontaneous curvature effects has not been pursued previously. Although the approximations in Eq. 32 do an excellent job for the lipid systems considered in this work, it is possible that other bilayers will exhibit stronger coupling between protrusion and bending modes. In such a case, it would be necessary to employ the full-forms (Eq. 26). In this work we use the expressions in Eq. 32 to fit simulation spectra and always find (see caption of Table 1) that the measure of coupling strength discussed above, 2²/(kk_c), is 0.3 or smaller. This provides us with a consistency check on the use of this approximation.

• They are continuous, rather than piecewise, expressions. This eases data fitting considerably, by avoiding a priori designation of individual data points to specific regimes. Although continuous expressions have been put forward (12) for the case of height fluctuations, no corresponding expressions have previously been discussed in the context of peristaltic fluctuations.
  
• The protrusion regimes of both spectra contain the constant term k. In the limiting case of q this constant is negligible, but, as shown in Appendix A, most large q-data for actual systems does not correspond to this limit.
  
• The monolayer spontaneous curvature contributes an effective negative surface tension to the thickness fluctuations (we find > 0 for all simulated data sets). This term, which was neglected in previous analyses, can result in a nonmonotonic thickness fluctuation spectrum.
  
• Due to our assumption of volume conservation, the compressibility modulus k_A appears in the thickness fluctuations (as opposed to a separate leaflet binding constant k_e (13–15)). In Comparison to Simulation Data, we discuss the validity of calling this constant k_A.
  
• We predict the height and thickness spectra to be characterized by the same bending rigidity. Furthermore, this bending rigidity corresponds to twice the bending rigidity of each leaflet. In analyzing their simulations, the authors of Lindahl and Edholm (13), Marink and Mark (14), and Chiu et al. (15) allowed thickness and undulation bending rigidities to differ from one another.
  
• The protrusion regimes of both the height and thickness fluctuations are identical. Although both Lindahl and Edholm (13) and Marrink and Mark (14) noted such equivalence in their data, they did not mandate it in their fitting procedure. In fact, according to the expressions in Eq. 32, the height spectrum should, in principle, contain no information beyond that already contained in the corresponding thickness spectrum. In practice, it is easier to extract k_c from the height spectrum than the thickness spectrum, so it is worthwhile to analyze both data sets simultaneously.
  

In summary, we have presented a theory with the same overall number of fit constants (five) as that presented in Lindahl and Edholm (13), but which includes the protrusion restoring force, the monolayer spontaneous curvature, and the area dependence of this spontaneous curvature.

	COMPARISON TO SIMULATION DATA

TOP ABSTRACT INTRODUCTION A GENERAL MODEL FLUCTUATION SPECTRA OF... COMPARISON TO SIMULATION DATA PROTEIN-INDUCED DEFORMATION... MOLECULAR SIMULATION DISCUSSION APPENDIX A: APPROXIMATIONS TO... APPENDIX B: CURVE-FITTING AND... ACKNOWLEDGEMENTS REFERENCES

 
Four sets of simulation data collected at vanishing tension were analyzed using the analytical model just described. The three data sets taken from atomic simulations (DPPC(13), GMO(14), and SM(15)) have appeared elsewhere. The coarse-grained data was taken using a model presented in Brannigan et al. (16) and readers are referred there for a full description of the model. The model parameters remain unchanged from that study, with the exception of the temperature. The simulations discussed in this article were all run at k_BT = 0.85, while simulations in Brannigan et al. (16) were run at k_BT = 0.9, where represents the energy scale = 2.75 kJ/mol. Systems at the lower temperature equilibrate more quickly. They also have a much lower monomer fraction (fraction of molecules that are not in the membrane), which simplifies analysis. Studies of phase behavior (not presented here) put the melting temperature of the five bead bilayer at k_BT = 0.7, so we are still well within the fluid region. Also, although Brannigan et al. (16) discussed a range of values for the molecular bending rigidity c_bend, all simulations in this article use c_bend = 7.

The height and thickness fluctuation spectra for DPPC, GMO, and SM are reproduced in Fig. 4, along with the corresponding data from our coarse-grained (CG) model. In the case of the atomistic studies, data presented in these figures necessarily reflect the conventions used in the original simulations. For DPPC, the position of the water-hydrophobic tail interfaces (z⁽¹⁾ + ⁽¹⁾ and z⁽²⁾ + ⁽²⁾ in our notation) were defined by the position of the carbon connecting the tails to the headgroup (13). A similar definition was used for SM (specifically, the location of the C13 atom (15)). In GMO, the interfaces were defined by the surfaces separating water from hydrophobic components of the bilayer core, irrespective of which carbon atom(s) along the chain happened to lie nearest the water (personal communication, S.-J. Marrink, 2005). In the CG model the interface was defined by the position of the second bead along the chain (the so-called interface bead (16)). In DPPC, GMO, and CG, the thickness is defined as one-half the distance between corresponding interface groups in opposing leaflets, whereas in SM the monolayer thickness was measured directly (the distance between leaflets and correlations between leaflets were neglected). It is worth emphasizing at this point that only the data collected for the CG model, GMO, and DPPC correspond well to the analytical model we have presented; SM data is somewhat different, due to the individualized treatment of monolayers.

For each system, we have two data sets (height and thickness spectra). We fit these data sets simultaneously to the expressions in Eq. 32 using the five elastic constants (k_c, k_A, , k, and ) as fit parameters. Our fitting algorithm is described in Appendix B and our fitting script is available on the Internet (38). Best-fit elastic constants are displayed in Table 1, along with the 95% confidence intervals for each fit parameter. The method for obtaining these confidence intervals is also described in Appendix B.

All four data sets are well described by the analytical model (expressions in Eq. 32), and, as demonstrated in Appendix A, the resulting numbers are consistent with the approximation that protrusion and bending modes are uncoupled. Furthermore, when we fit these data sets to the full form for the spectra (expressions in Eq. 26), we obtain curves that are essentially identical to those resulting from the approximate fits, and the fit parameters all agree within confidence intervals. As also shown in Appendix A, the effective negative tension induced by finite positive can play a significant role in determining these spectra. In fact, this negative tension provides a possible explanation for the nonmonotonic behavior (previously (13) attributed to poor sampling) seen in two of the atomistic spectra.

The bending rigidities for undulations remain essentially unchanged from the original analyses, but we report different values for the remaining constants. Reported values for the interfacial tension originally ranged from 8.5 mN/m² for GMO (14) to 40 mN/m² for SM (15). Since this number is expected to primarily represent the tension associated with the oil-water interface (with some renormalization due to lipid shape (11)), such variation is surprising. Upon reanalysis, we report ranging from 13 mN/m for GMO to 23 mN/m for SM, which is a more reasonable range. Upon correcting for molecular shape as advocated by Lipowsky and Grotehans (11), we obtain a range of 11–18 mN/m. Theoretical predictions suggest that should fall between 20 and 50 mN/m (39).

The compressibility modulus k_A is usually measured by area dependence of the surface tension (12) or fluctuations in the area per molecule (40). The latter method was used to measure k_A in the three studies already published (DPPC, GMO, SM) and, consequently, we have also used it to measure k_A in the coarse-grained model. Table 1 compares k_A as measured by the area fluctuations, to k_A as measured by the thickness spectrum, for all four systems studied. There is clear order-of-magnitude correspondence, and the discrepancies do not appear systematic: for DPPC and SM the area fluctuations suggest a greater k_A, whereas for GMO and CG, the thickness spectrum suggest a greater k_A. (We assume it is a coincidence that DPPC and SM are the double-chained lipids, although GMO and CG are single-chained.) Aside from finite size-effects and noted problems in converging k_A (40), other possible sources of the discrepancy are that:

1. The volume per lipid is not strictly conserved in real systems. In this case, our analytical model is still valid, but it is not appropriate to identify the constant we have been referring to as k_A with the true area compressibility modulus.
   
2. Our assumption that protrusions do not enter into the volume conservation condition may be incorrect. If the membrane can pull volume from as well as z, one expects our measurement of k_A via the thickness fluctuations to be too high. Protrusions allow for another degree of freedom to facilitate lipid compressibility. Since we do not see a systematic trend of thickness determined k_A being too large, this possibility cannot account for all observed discrepancies.
   
3. Definition of the bilayer thickness is subject to a certain amount of ambiguity as evidenced by the various procedures employed in calculating this field by the original workers. For instance, in GMO, the thickness spectrum reflects a more general definition of the interfacial surface than was introduced in the other simulations. Although the interfaces defined for GMO probably correspond more precisely with the true interfacial surfaces than the other simulations, this definition actually contrasts somewhat with our theoretical treatment of area/lipid, conservation of hydrophobic volume, and related issues. It is possible that the procedure used for collection of GMO data tends to underemphasize thickness fluctuations, relative to the other studies. Such an effect would lead to an anomalously high value for k_A. Similarly, in the SM data, each leaflet has been Fourier-transformed and squared separately; consequently, this data neglects correlations between the leaflets. If thickness fluctuations are anticorrelated across the two leaflets, then this data would result in a measurement for k_A that is too low. (These variations in collection methods of the original MD data could also potentially affect the values obtained for other elastic constants. Since k_A is the constant most directly responsible for thickness deformations, it is reasonable that we see the largest discrepancies here.)
   

As previously stated, we cannot directly extract c₀ or c₀' from the thickness fluctuations; we can only measure , which depends on both. Since fluctuation analysis does provide the value of monolayer bending rigidity, k_c/2, we may infer c₀ from analysis of the stress profile (1),

(33)

where (z) is the surface tension density at height z and the integration is over the whole bilayer, centered at z = 0. The stress profile (Fig. 3) was measured for the CG model, as documented elsewhere (16,41,42). To avoid the smoothing effect of undulations (16), we measured the stress profile in a small system (128 molecules). Using this method (and k_c obtained from the spectra), we estimate c₀ = 0.041 nm^–1. Using the fluctuation spectra values /t₀ = 0.1 nm^–2, t₀ = 2.4 nm, and ₀ = 0.59 nm², we can further estimate c₀' = –0.34 nm^–3. This information will be useful in comparing to the information obtained from the protein-induced deformation data of the following section.

	PROTEIN-INDUCED DEFORMATION PROFILES

TOP ABSTRACT INTRODUCTION A GENERAL MODEL FLUCTUATION SPECTRA OF... COMPARISON TO SIMULATION DATA PROTEIN-INDUCED DEFORMATION... MOLECULAR SIMULATION DISCUSSION APPENDIX A: APPROXIMATIONS TO... APPENDIX B: CURVE-FITTING AND... ACKNOWLEDGEMENTS REFERENCES

 
Theory
Consider the bilayer model discussed in Fluctuation Spectra of Homogeneous Membranes, above, but with a rigid cylindrical inclusion of radius R and thickness 2D embedded in its center. More precisely, 2D represents the thickness adopted by the membrane right at the edge of the protein. This value is presumably set by the presence of hydrophobic and interface favoring residues around the protein's exposed surface. We seek to predict the deformation profile t(r), where r is the distance to the inclusion center (Fig. 5). Since height and thickness deformations are uncoupled, we need only consider the second and third lines of Eq. 21:

(34)

Instead of integrating over the whole sheet, as in the previous section, our domain of integration is a square of side length L with a circle of radius R cut out of the center to accommodate the protein. Assuming that the inclusion cannot tilt, the membrane surface is subject to the boundary condition at r = R of

(35)

It should be recalled that t₀ is the equilibrium monolayer thickness, whereas t is defined as the deviation in monolayer thickness away from this value.

We imagine ^(1,2) to represent microscopic noise or roughness at the lipid interface. With this interpretation, it is reasonable to assume that our thickness condition above translates to ^–(R) = 0 and z^–(R) = D – t₀, in terms of our field variables. Because protrusion and bending deformations are coupled, this condition is not equivalent to setting ^–(r) = 0 for all r, even for purposes of computing the thermally averaged deformation. As described below, however, the inclusion deformation profile for our CG bilayer is consistent with either constraint ^–(R) = 0 or ^–(r) = 0 everywhere since protrusions and bending modes are very weakly coupled in this system. If (r) is set to 0 for all r, Eq. 34 is equivalent to the free energy minimized analytically by Aranda-Espinoza et al. (31) to predict the zero temperature membrane deformation profile for a two-dimensional array of inclusions.

We solve for the average deformation profile predicted by Eq. 34 using Metropolis Monte Carlo simulation on a discretized periodic square lattice. Each lattice site has an associated value of z^– and ^–. A trial move consists of randomly choosing a lattice site, which field to perturb, and the perturbation size. The lattice has either 40 x 40 sites with a spacing of 0.75 nm (CG parameters) or 60 x 60 sites with a spacing of 0.25 nm (GMO parameters). These choices for the number of sites and spacing were verified to be sufficiently large and small, respectively, that further refinement and system expansion did not alter our results. Lattice sites within a radius R from the lattice center were discarded, since the membrane does not exist inside the inclusion radius. The only conditions directly imposed during the Monte Carlo were periodic boundary conditions and ^– = 0 and z^– = t_R at the inclusion-lipid boundary. This means that membrane slope and curvature at the inclusion/lipid boundary are unrestricted and protrusions are allowed to occur everywhere except at the sites immediately adjacent to the inclusion.

We ran the lattice simulation for two sets of elastic parameters: those resulting from the CG simulations and those resulting from the GMO simulations (Table 1). These correspond to a system in which ^– and z^– are very weakly coupled and a system in which ^– and z^– are more strongly coupled, respectively. For the spontaneous curvature, which cannot be obtained from the fluctuation spectrum, we either used a value obtained from the stress profile (CG) or assumed that c₀ = 0 (GMO). In both cases we set R = 2.25 nm and t_R/t₀ = 0.094, which coincide with the radius and the height, respectively, of the nontilting inclusion used to induce a deformation in the CG membrane (described in Molecular Simulation, below).

The numerical results for CG are shown in Fig. 6 A. Although ^– does not vanish for all r, neither does it ever constitute an appreciable amount of the deformation. Both z^–(r) and t(r) = z^–(r) + ^–(r) very nearly obey the analytical expression derived in Aranda-Espinoza et al. (31). Analytical treatments based on Hamiltonians similar to Eq. 34 do not directly solve for the thermally averaged deformation profile (22,29–31,43–46). Instead, the minimum energy configuration is determined by taking the functional derivative of F_inc (in our example), with respect to z^–(r), assuming ^– = 0 everywhere and cylindrical symmetry of the solution around the inclusion (31). The resulting expressions for the minimum energy deformation profile are fourth-order differential equations, the order insured by the bending (k_c) terms in the Hamiltonian. In our Monte Carlo treatment, only two boundary conditions are necessary to run the simulation as discussed above and one of these relates to the protrusion modes that are neglected in analytical treatments. For the purposes of minimization, however, we need a total of four boundary conditions to solve the differential equation. In practice, two of these boundary conditions typically relate to the deformation profile at points distant from the inclusion (to keep the solutions bounded at infinity (22), or to enforce symmetry constraints within periodic boundary conditions (31)) with two boundary conditions to be specified at the edge of the inclusion. One obvious boundary condition simply sets t(R) = D as we have introduced above.

View larger version (9K): [in this window] [in a new window]   FIGURE 6  Average deformation profiles of the fields z– (solid) and