|Home | About | Journals | Submit | Contact Us | Français|
One of the many aspects of membrane biophysics dealt with in this Faraday Discussion regards the material moduli that describe energies at a supramolecular level. This introductory lecture first critically reviews differences in reported numerical values of the bending modulus KC, which is a central property for the biologically important flexibility of membranes. It is speculated that there may be a reason that the shape analysis method tends to give larger values of KC than the micromechanical manipulation method or the more recent X-ray method that agree very well with each other. Another theme of membrane biophysics is the use of simulations to provide exquisite detail of structures and processes. This lecture critically reviews the application of atomic level simulations to the quantitative structure of simple single component lipid bilayers and diagnostics are introduced to evaluate simulations. Another theme of this Faraday Discussion was lateral heterogeneity in biomembranes with many different lipids. Coarse grained simulations and analytical theories promise to synergistically enhance experimental studies when their interaction parameters are tuned to agree with experimental data, such as the slopes of experimental tie lines in ternary phase diagrams. Finally, attention is called to contributions that add relevant biological molecules to bilayers and to contributions that study the exciting shape changes and different non-bilayer structures with different lipids.
I fondly recall Faraday Discussion 81 in 1986 on lipid vesicles and membranes, and I am pleased that the Faraday Society has chosen to revisit this vitally important area. That early Discussion brought out some issues that needed to be frankly discussed. The Faraday Discussions are especially valuable because such exchanges are recorded for those not able to be present and who would otherwise not be privy to candid points of view. I hope that the present Discussion will also not shy away from the controversial and I will do my best, in this introductory lecture, to initiate such discussion. Similarly, Barry Ninham in his introductory lecture to Discussion 81, while praising progress, compared the then present state of the field to the satirical novella Flatland, and he then proceeded to elaborate on his list of unresolved problems.1 At the time, I thought Barry might have overdone it, and some of you may think that of my lecture for this Discussion. On the other hand, in a wise and influential book, Efraim Racker entitled his first chapter “Troubles are good for you”.2 He also attributed to Clerk Maxwell the quote “if you have built a perfect demonstration do not remove all traces of the scaffolding by which you have raised it.” Especially in view of Maxwell’s equations building on Faraday’s scaffolding, the Faraday Discussions format appropriately provides an opportunity to examine scaffolding.
However, before engaging with specific troubles, it is important to emphasize for those outside the field that there is much solid chemical and physical information about biomembranes and there are valuable theoretical paradigms that have focused on relevant basic quantities. We know that the most important biomembranes are fairly flexible structures only a few nanometers thick with a basic core of lipids containing a large variety of biochemically interesting proteins. While this general knowledge already suffices for much of general biology, deeper understanding of why there are such significant compositional differences between the many different membranes that occur in biology and even among the different membranes in a typical cell is required and physical studies must extend from the molecular to the supramolecular level. Without further ado, let me dive into what I think are a few of the outstanding fundamental problems, in the hope that this will promote further advancement of the field.
It is well recognized, as evidenced in papers at this Discussion,3–8 that the bending modulus KC is a most important quantity for describing the flexibility of membranes at a mesoscopic, supramolecular level. Other material quantities are the Gaussian curvature modulus KG,7,8 important for membrane fusion,9 and the isothermal area compressibility modulus KA.10 Nevertheless, the centrality of the bending modulus is indicated, in the case of the Gaussian curvature, by the fact that it is the ratio KG/KC that is usually reported rather than KG alone, and the compressibility modulus KA is also often closely related to KC (vide infra). Also, there are relatively few measurements of KA and KG, so let’s focus on KC, without in any way disparaging the importance of KG and KA that are central to several papers in this Discussion.
For a symmetric membrane, the local (free) energy of bending is universally approximated as
where C(r) is twice the mean curvature at location r = (x, y) in the plane of the membrane that has fluctuations in the z direction. Let me emphasize that differences of a factor of two in KC can be highly significant, especially when considering thermally activated processes because rate constants depend exponentially on the activation energy. For example, for lipid membrane fusion, activation energies for curved fusion structures greater than 40kT are deemed to be kinetically incompetent at the biological time scale,11,12 so it is important to make estimates of the energies of intermediates predicted by any theory.13 Such an estimate might be 30kT if one measured value of KC is used. However, if another measured value of KC is twice as large, then the calculated activation energy would double, the predicted rate constant would decrease by 13 orders of magnitude, and the theory would be kinetically incompetent. Even if theories rely on energetically driven processes rather than thermal activation, factors of two in KC could still discriminate whether a theory is energetically competent or not.
The trouble is that reported values of KC do indeed differ by factors of about 2. Table I shows some of the notable more recent values of KC/kT for three of the most popular lipid bilayers, DOPC, SOPC and DMPC, as well as for the thicker diC22:1PC bilayer. In some cases, values for a particular lipid differ considerably for the same method of measurement. Furthermore, there are differences when comparing values obtained from different methods of measurement. The most common methods shown in Table I are abbreviated: (1) diffuse scattering (X-ray) on stacks of membranes (X-ray stacks), (2) micromechanical manipulation (MM) of giant unilamellar vesicles (GUVs), often called the aspiration pipette method (MM GUV), (3) pulling cylindrical tethers from GUVs (tethers GUV) and (4) analysis of shapes of GUVs (shapes GUV). In a review, Derek Marsh separated values of KC into three tables according to methodology, in order to focus on differences between different lipids and not dwell on the differences between methods.14 I take the opposite perspective in Table I in order to focus on the fundamental issues of KC measurements.
First, note that the temperatures of the measurements cited in Table I often differ. Based on X-ray data for DOPC,15 the factor for the temperature dependence of KC is only about 0.995/°C and this factor has been used to adjust all the literature values to T = 30 °C. Next, regarding the X-ray column, my group has taken repeated data on these lipids and the values shown here are updated from several papers and from unpublished results. We often quoted uncertainties of ~8% in papers, but when subsequent data sets have been acquired and analyzed by different co-workers, the uncertainties should probably be raised closer to 15%. Two other groups have also analyzed diffuse X-ray scattering in different ways to us. Agreement for DOPC is satisfactory16 and it was satisfactory for DMPC17 until the erratum published by one of the groups.18
I have found it satisfying that our X-ray results agree so well with the earlier MM results from Evan Evans’ group.19 However, agreement is not so good among the different groups that have subsequently used MM. One source of disagreement stems from the sugar solutions employed to improve visual contrast. Although a priori I would have supposed that the 200 mM sucrose/glucose employed in ref. 19 would be innocuous, two groups have reported a strong dependence on the concentration cS of sugar.23,20 To emphasize this, for these two references I have entered two values of KC/kT in Table I for SOPC23 and DOPC.20 The larger value is for no sugar and the smaller value is for cS = 200 mM; these values were obtained from the original papers as follows: Vitkova et al. combined values of KC for SOPC with cS in the range 100–300 mM using the MM method with values obtained using the shape method with cS = 0 and 50 mM from which they produced a formula, KC(cS) = KC(0)exp(−cS/200 mM).23 Using their fitting formula gives KC/kT = 30 for a zero sugar concentration and KC/kT = 11 for 200 mM sugar, which are the two numbers shown in Table 1. For DOPC, a decrease in KC/kT from 23 in 8 mM sugar to KC/kT = 12 in 100 mM sugar was reported using the MM method;20 assuming that the decrease with increasing cS is also exponential, extrapolating to 200 mM sugar gives the value KC/kT = 6 that I show this in Table 1. This last value for DOPC and the value KC/kT = 11 given in the table for SOPC in 200 mM sugar are much smaller than the values reported for cS = 200 mM by Rawicz et al.,19 and the corresponding values given for cS = 0 are larger.
Another source of disagreement involves a different method of analysis of MM data that employed global fitting of the MM data to both KC and KA for all tensions,24 instead of separately analyzing the low and the high tension regimes.19 This gave the largest value of KC/kT = 36.9 for SOPC in the table.24 Interestingly, 75 mM sugar was used.24 Employing the exponential formula above23 then would predict KC/kT = 19.8 for 200 mM sugar, in good agreement with Rawicz et al.19 However, the formula also predicts KC/kT = 51 with no sugar, which is far higher than the value reported by Vitkova et al.23 These effects of sugar concentration and MM analysis methodology are intriguing. It would be very interesting to obtain KC/kT for increasing sugar concentrations using one of the other methods. In the meantime, I will take the values from the Evans’ group19 as the canonical MM results. One rationale for this choice is that their values are beautifully understood in terms of their polymer brush model, which gives the number 24 in the following formula relating KC to the area compressibility modulus, KA
where 2DC is the thickness of the hydrocarbon core of the bilayer.19 My group also subjected the polymer brush theory to additional tests involving the measured temperature dependence of the area/lipid A and it passed our tests with flying colors.15
Shape analysis of giant unilamellar vesicles (GUV) is the oldest method in Table 1. As this method involves direct visualization of fluctuating bilayers, it would seem to be unequivocal. However, early values varied widely and considerable efforts have been made to improve the method.34,27,33 While the results in Table 1 for SOPC indicate some differences using this method, the broader conclusion often noted in the literature14,22,24,33 is that values for KC/kT obtained from the shape method tend to be larger than those obtained from MM. The last column in the table roughly quantifies the ratio of the values from these two methods. The relatively few values that I found from the tether method21,22,25 do not exhibit clear trends relative to the MM and shape values.
I would like to speculate on why values of KC might legitimately be larger from shape measurements than from MM measurements. This speculation begins by noting the different length scales of the experiments. The length scale L for shape analysis is typically of the order 10μm (wavelength of the 12th mode in a GUV with radius R =20 μm). In contrast, what matters in MM is the pushing out by aspiration pressure of the excess area due to undulations over all modes with wavelengths down to the order of a = 0.8 nm, the distance between lipid molecules. From eqn (1), the tensionless excess area is derived to be proportional to (1/KC)log(R/a), where R/a spans 4 or more decades of wavelength. Because of the logarithmic dependence, each decade contributes equally to the excess area, while only the top decade contributes to the shape analysis because smaller wavelengths are not optically visible. This distinction in the length scales makes no difference in extracting KC if the energy is given by eqn (1). That energy gives a height–height undulation spectrum with intensity I(q) ~ 1/(KCq4) varying inversely with the 4th power of wavenumber q.
However, there is an alternative theory35 that adds an energy that depends on molecular tilt t(r),
It should be noted that a term inversely proportional to q2 was ascribed to protrusions for over ten years, but it has been emphasized that there is no theory of protrusions that gives the same functional form as eqn (4).36 Also, recent analysis of a simulation concluded that any contribution of protrusions to I(q) was small compared to the effect of tilt.37
If eqn (4) is true, then the q−4 term would dominate I(q) at the large length scale and KC would be obtained from shape measurements. The added term due to tilt in eqn (4) would become important as q2 exceeds Kθ/KC and this would add an extra contribution to the area pulled out by aspiration-induced tension. When analyzed assuming eqn (1), the apparent KC obtained from MM would then likely be smaller than the KC in eqn (4). Such a reduction in KC is qualitatively consistent with the analysis of simulations;36,38 reduction factors ranging from 1.4 to 4 depend upon the simulation and how it is analyzed. Indeed, the method of simulation analysis I have preferred shows no q−2 term at all.39,40 That would be consistent with tilt simply not involving the energy term in eqn (3) for bilayers in the chain melted, fluid phase. However, molecular tilt does exist and a chain tilt order parameter is even obtained by a different kind of X-ray method41 than the one used to obtain KC/kT.
What about the KC values from the X-ray method? These were obtained from stacks of bilayers that have an interbilayer interaction energy that effectively constrains the in-plane length scale to the range of wavelengths smaller than 0.5 μm. As the analysis also assumed only eqn (1) for the individual bilayers, this speculation suggests that the X-ray values for KC/kT in Table 1 should indeed be smaller than the values obtained from the shape method. It is a bit surprising that the X-ray values are close to the MM values as each would presumably have a different, as yet unknown, reduction factor. Indeed, I have not yet given up on the possibility that the X-ray and MM method give the true KC values and that the shape values are artifactually large.
Nevertheless, if it is decided that this speculation is valid, then the obvious way forward is to analyze MM and X-ray data by adding a tilt energy. Of course, this adds another parameter Kθ that may not be reliably extractable from either MM or X-ray data. It would also involve heavy theoretical development for the X-ray method, which, like the other methods, is already quite complicated.42 In the short term, it may be more practical to continue to report apparent KC values just using eqn (1) in the MM and X-ray analyses, recognizing that these apply to smaller length scales than those that are visible in the shapes of GUVs. Indeed, it may also be argued that it is better for many uses in biophysics to limit the set of descriptors and not proliferate parameters by including Kθ. In that case, if one wishes to evaluate the kinetic competence of thermally activated theoretical intermediates using KC from MM and X-ray measurements, then one should not invoke a tilt degree of freedom. If one does invoke tilt,11 then the value of KC should be increased to the larger values determined by the shape measurements.
Before finishing with the KC story, let me mention a very interesting effect of cholesterol on the stiffness (i.e., KC) of membranes. It had been the conventional wisdom that cholesterol stiffens fluid phase bilayers. Indeed, it does stiffen DMPC, which has two saturated hydrocarbon chains, but it doesn’t stiffen DOPC, with two mono-unsaturated chains, at all.30,43 Cholesterol stiffens SOPC30,43 and POPC,44 both with one unsaturated chain, but less than half as much as it stiffens DMPC. The remarkable result that cholesterol does not stiffen DOPC bilayers has been amply confirmed.33,21,22 Even though different investigators using different methods report different KC values for pure DOPC, the null cholesterol dependence in DOPC is found using several methods, so this appears to be a robust result. This warrants asking, what are the physical chemical interactions responsible for this striking difference in the effect of cholesterol on different lipid bilayers? It has been reported that cholesterol interacts more strongly with DPPC than with DOPC,45 and it would be interesting to develop a theory that would use that result to predict the KC effect quantitatively.
Even if we don’t understand the physical chemistry of the cholesterol effect very well, it is interesting to connect the effect of cholesterol on the KC value for different lipid bilayers to the well known biomedical assessment that saturated fats are not as healthy as unsaturated fats. Despite its pernicious effect in the bloodstream, cholesterol is required in mammalian cells for many purposes.46 However, cell membranes are required to be flexible; while stiffer membranes could be viable, they would incur slower activation rates and impose larger metabolic loads in order for cell membranes to undergo processes that require curved intermediates, such as endocytosis and exocytosis. Both these desiderata, flexibility and the presence of cholesterol, are provided by DOPC bilayers composed of unsaturated fatty acid hydrocarbon chains, but not by DMPC bilayers composed of saturated fatty acid chains.
Finally, the goal of having some simple relation between KC and KA, like the one in eqn (2), is desirable in that it reduces the number of independent quantities necessary to describe membranes. This goal is strongly supported by results for many pure bilayers for which the polymer brush theory works well,19 but it breaks down when cholesterol is added.43 This is not surprising as the rigid cholesterol ring structure should modify the polymer brush model. However, this emphasizes that KC and KA really are two independent quantities whose relation depends on molecular details beyond the thickness of the membranes. Clearly, the cholesterol story remains as fascinating as ever. Not surprisingly, cholesterol appears in many papers in this Discussion.47–54
This subsection serves as a connector between KC and the next section that will focus on membrane structure. As was shown earlier,19 our data in Fig. 1 show that KC is well correlated with the structural thickness DHH. The temperature dependence of DOPC is also shown, and that was used to adjust the DPPC thickness. According to simple elasticity theory, the KC dependence should be quadratic relative to some core thickness, often assumed to be the hydrocarbon thickness, if all other factors are equal. Despite the difference in unsaturation of the hydrocarbon chains, a quadratic formula does fit well for these PC lipids. The best fit is to (DHH − 7.4)2 and the offset of 7.4 Å suggests that the boundary of this core is located 7.4/2 = 3.7 Å closer to the center of the bilayer than the phosphate group, which is located close to DHH/2 (but note the caveats in the next section).
It is universally understood by biophysicists that one needs structure in order to best understand function and I will not belabor this important perspective. It is also recognized by lipid biophysicists that the lipid content matters for function, and this motivates the goal of obtaining enough structural accuracy to determine differences between bilayers composed of different lipids. As is well known, a faux trouble with obtaining membrane structure is simply that most biomembranes exist in the hydrocarbon chain melted, dynamically flexible, fluid states, so the kind of structure that small molecule crystallography obtains at sub-Angström level resolution just isn’t there. This is conspicuously obvious when one examines snapshots of simulations. Statistical averaging over many snapshots and many lipids, as in Fig. 2, shows that the phosphate groups are distributed in the direction along the bilayer normal with full width at half maximum of the order 6 Å. Nevertheless, this really is a Definite structure, but at the sub-nano scale instead of at the sub-Å scale. Indeed, the electron dense phosphates mainly (but not quite, vide infra) determine the so-called head–head thickness DHH that I used in Fig. 1. DHH is the best determined structural parameter that X-ray methods provide. X-ray methods also provide the repeat spacing D for multilamellar samples, but D only gives the sum of the water thickness DW and the volumetric bilayer thickness DB; disentangling these using a Luzzati-type analysis has been fraught with difficulties.56 A much better way to obtain the bilayer thickness DB employs neutron scattering on unilamellar vesicles.57,58
The basic X-ray and neutron data are the bilayer form factors F(qz) that are the Fourier transforms of the electron density and neutron scattering length density, respectively. I will leave experimental and analytical details to primary papers42,60,61,32,57 and a forthcoming review that is under preparation. It suffices here only to emphasize that the X-ray data in Fig. 3 were obtained from fully hydrated oriented multilayered samples with ample water space between the bilayers (these are the samples that the KC values in the X-ray column of Table I were obtained from) and also from unilamellar vesicles (ULV); agreement is good between the two kinds of samples for qz values where the data overlap.32,31 The neutron data were taken exclusively on ULV.57,58 All these data are continuous in qz, unlike the discrete orders obtained from multilamellar samples at reduced humidity when the headgroups are jammed into contact with no clear water cushion.
While the X-ray data provide DHH and the neutron data provide DB rather directly, one wishes to know the distributions of different parts of the lipid molecule along the bilayer normal, as well as the distribution of additives, such as cholesterol and peptides. The traditional route to these quantities employs modeling that incorporates as much other information as possible. In particular, measured volumes56 provide a connection between distributions along the bilayer normal and the in-plane area of lipid62 (plus additives if present43,63,64); this area A is the main structural information along the surface of a membrane for single component lipid bilayers. The results of modeling provide reliable electron density profiles and neutron scattering length density profiles. The trouble is that there are many constraints that have to be employed in modeling59,65 in order to obtain the more detailed component distributions shown in Fig. 2 that one would like to see.
Another route to obtaining the detailed distributions of components is to use simulations. Of course, even assuming that the simulation is competently carried out and that there is enough computer time to equilibrate the simulation,66 there is still trouble regarding whether the force fields are valid. An early, apparent success was the outstanding agreement for the electron density profile and the X-ray F(qz) for DOPC using CHARMM27.67 The only snag at that time was that, when run at zero surface tension, C27 bilayers are too thick and have too small area/lipid A; consequently, their F(qz) values do not agree with experiment. However, I continue to believe that this area shrinkage could be due to the water model being too primitive and so using a non-zero surface tension or a constant area simulation is justified,65,68 contrary to the often stated goal of obtaining agreement with experiment using tensionless NPT simulations.69,70 The outstanding agreement with X-ray F(qz) seen in Fig. 3 was obtained when the area was fixed to the value A = 72.4 Å2 that agreed with the value obtained from modeling, so the simulation and the model are mutually supporting.
The wonderful agreement in the preceding paragraph has been destroyed by adding the requirement that the simulation must also agree with neutron data, Fig. 4. This first became apparent, although it was not emphasized, in a paper that focused on modeling and not on testing simulations.59 The goal was to develop modeling that could use both X-ray and neutron data simultaneously to provide enhanced structure. Indeed, that development used the C27 simulation of DOPC with A = 72.4 Å2 to inform the model building and to validate the final model that was developed. The first trouble came from the values of A obtained by applying the model to real data. Although A for DPPC did not change much from earlier modeling based only on X-ray data,55 the area for DOPC dropped from 72.4 Å2 to 67.5 Å2. The second trouble, which was not emphasized,59 is that CHARMM27 did not agree with both the neutron and X-ray data for DOPC. When run at A = 72.4 Å2, it agreed with the X-ray data but it disagreed with the neutron data, as shown by the red lines in Fig. 3 and and4.4. When run at A = 67.5 Å2, C27 agreed with the neutron data, but disagreed with the X-ray data, as shown by the blue lines in Fig. 3 and and4.4. (The data and the simulation can also be viewed more closely from the website http://www.norbbi.com/public/public.html#Software using a software package SIM-toEXP especially designed for comparing simulations with scattering data.71)
The first trouble with obtaining A using only X-ray data has been located. Although the DHH thickness is obtained rather directly from X-ray data, to obtain A we had divided the hydrocarbon volume VC of the lipid by DC, where 2DC is the hydrocarbon thickness, but DC had to be estimated. The difference (DHH − 2DC)/2 defined as DH1 was found to be about 5.0 Å, both from the DMPC gel phase structure72 and from C27 simulations, so this value of DH1 was assumed to be the same for all lipids.56 However, modeling to both neutron and X-ray fluid phase data resulted in a large decrease to DH1 = 3.7 Å for DOPC; although, for DPPC, a much smaller decrease to DH1 = 4.7 Å was obtained. While there are far fewer neutron data than X-ray data, only extending to about 0.2 Å−1 in q, as seen in Fig. 4, compared to 0.8 Å−1 for the X-ray data in Fig. 3, the neutron data suffice to determine from DB 61 which A is directly obtained by dividing DB into the lipid (or lipid plus additive) volume VL. Clearly, future determination of A should include both neutron and X-ray data.
Although the trouble involving CHARMM simulations has been known since 2008, it is not yet resolved. Let me give a progress report here with some speculation about what the trouble may be due to. First, note that the current version of CHARMM is C36, which is improved over C27 in that it obtains better values for the NMR order parameters near the top of the hydrocarbon chains and in the glycerol group.73 It should, of course, be reiterated that NMR order parameters are crucial data that have been used for a long time to test simulations. Agreement with the NMR data is valuable primarily to test the hydrocarbon chain force fields, particularly when the more modern NMR data come from perdeuterated chains, which do not allow identification of the carbon position to the order parameters at the top of the chains.74 C36 also obtains areas in tensionless NPT simulations that are closer to the range given by experiment. Unfortunately, C36 is slightly worse than C27 for obtaining agreement at the same A with both the DOPC raw neutron and X-ray F(qz) data (unpublished), although it agrees reasonably well for saturated DPPC, DMPC and DLPC.73
Fig. 5 schematizes the general outcome that will occur for any simulation method that is run at many different areas in NPAT simulations or equivalently for many surface tensions γ in NPγT simulations. The simulation will best fit experimental neutron data at Aneutron; it will best fit X-ray data at AX-ray, and it may best fit other data, such as NMR order parameters, at some Aother. Of course, the ideal simulation should have
I believe that satisfying eqn (5) is a test that should be applied to any force field and for a variety of lipids. Furthermore, I will next argue that the information obtained by applying this test may be instrumental in improving the force fields.
In contrast to most NMR data that focus on the hydrocarbon chain region, neutron and X-ray data focus on the interfacial, headgroup region. Combining both sets of scattering data with simulations promises to elucidate this complex region. Although my focus until recently was on the DH1 difference between hydrocarbon thickness and the peaks in electron density profiles, I now suggest that a more revealing quantity is the difference between the volumetric thickness DB and the X-ray thickness DHH. Results for this quantity are shown in Fig. 6 as a function of A from experiment and from various simulations. (I thank the many simulators who have put their simulations into the SIMtoEXP format71 from which I then obtained the results shown.) The figure includes four popular lipids, DMPC, DLPC, DPPC and DOPC. Two different areas are shown for the experimental results for each lipid. The newer modeling using both X-ray and neutron data gives smaller areas compared to the older areas obtained from X-ray data only. However, for discussing the experimental DB − DHH, the differences in experimental areas are inconsequential because DB − DHH is essentially constant as a function of A.
In order for a simulation to agree with both neutron and X-ray data, it has to agree with both DB and DHH. By adjusting A, a simulation can be made to agree with either DB or with DHH, but to agree with both at the same A, it also has to agree with the difference DB −DHH. Fig. 6 shows that DB − DHH steadily diverges from the experimental values as A increases for C36 simulations, regardless of whether one compares different lipids or the same lipid at different areas (DPPC or DOPC in Fig. 6). This is consistent with the result that C36 agrees better with both X-ray and neutron F(qz) data for DPPC than for DOPC. In contrast, GROMACS (GMX) united atom (UA) simulations kindly run for me by Anthony Braun (DOPC and DMPC) and some others by Olle Edholm (DMPC) give similar values for DB − DHH as the experiments, as also shown in Fig. 6. Consistent with the above prediction, the GMX simulations give good agreement with both the DOPC X-ray and neutron data at the same area A; furthermore that area is A = 67.5 Å2, which agrees well with modeling.59 In an earlier paper that had not compared to neutron data, an NPT GMX simulation also gave good agreement with DOPC X-ray data for A = 67.5 Å2.75 GMX also gives good agreement for both neutron and X-ray data for DMPC at A = 60.8 Å2; although the agreement with DMPC data is slightly better for the all atom C36 simulation. It may be noted that using the Berger force fields76 in GMX does have one serious flaw, namely, the ratio r of the volume of the terminal methyl to the methylenes on the hydrocarbon chains is larger (r ~ 2.7) than experiment (r ~ 2),62,77 resulting in the methyl trough in the electron density profile being too deep. Indeed, the GMX simulation of DMPC with A = 62.2 Å2 in Fig. 6 used a different force field that gives r values close to 2.78 However, the r value has a relatively small effect on the F(qz) in the fluid phase.
At this point one might think that CHARMM simulations should be ignored in favor of GMX simulations. However, with even more detailed comparison of the different simulations and experiment, it should be possible to determine better what works and why. Note that the trouble may not just be the force fields. A recent paper reports that conformational equilibration in the headgroup region is exceedingly slow, even for single lipid molecules not in a bilayer.79 It would be good to learn the effect of the initial state of the headgroup conformations on simulated F(qz) and DB − DHH.
A plausible working hypothesis is that differences between GMX and C36 simulation results are due to conformational differences in the headgroup region. I can not address that question because simulation data in the SIMtoEXP format71 do not allow analysis of molecular conformations, but they do allow some other observations. Even though the location of the electron density profile peak at DHH/2 is usually identified with the average position zP of the phosphate, DHH/2 may be shifted to positions smaller than zP by the carbonyl/glycerol (CG) moiety. The shift in zP − DHH/2 will be larger if the force fields give CG a smaller volume that results in a larger electron density. The SIMtoEXP software has an app that obtains component volumes from simulations.80 It shows that the CG volume is indeed smaller for GMX and that, for both DOPC and DMPC, zP − DHH/2 is ~0.6 Å larger for GMX than for C36. This would account for the differences in DB − DHH at low A in Fig. 6 without invoking conformational differences. It may also be noted that for DOPC the shift zP − DHH/2 is close to zero for C36 but increases to ~0.3 Å for DMPC and DPPC. Of course, the location of the maximum of the total electron density also depends upon the interfacial water, which is usually modeled by SPC in GMX and TIP3P in CHARMM. Another suggestion is that the new, smaller values of DH1 that are required in the modeling of DOPC data may be due to greater tilting in the carbonyl–glycerol backbone region so that the phosphate distribution would lie closer to the Gibbs dividing surface of the hydrocarbon region even without molecular conformational differences. However, the quantity zP − DC is not much different for GMX versus C36 and for DMPC versus DOPC; so tilt and conformational differences are not presently indicated by simulation data in the SIMtoEXP format, but these should, in future, be extracted from the full simulation trajectories.
Despite these troubles that have persisted for some time now, I am optimistic that even better quantitative structures of simple lipid bilayers will be obtained in the not too distant future by close interaction between experimentalists and simulators. Here, I would like to address how I think this interaction should proceed. A decade ago, I was confident that we were able to obtain area A from X-ray scattering.56 Accordingly, area was a focal point of the interaction between simulation and experiment. However, because of the crudity of the non-quantum mechanical water potentials, it seemed better to use the experimental values of A to guide the simulation, either using NPAT or NPγT. Now that we are less confident in the values of the areas, both from scattering experiments and, even earlier, from NMR order parameters,81 it seems clear that, while area A remains a most important focal quantity, the development of simulations should not regard the goal of matching uncertain experimental A values as a primary goal.† Instead, simulations should compare to raw experimental data that is not subject to the modeling that is required to obtain A.63,65,75 If that comparison is satisfactory, then a model free estimate of A can be obtained from the simulation by finding the value of A that best fits the experimental data.65 In this approach, the surface tension γ is a one parameter work-around to compensate for the poor water–lipid interaction potentials.
Although I have focused here on raw X-ray and neutron F(qz) data, we also obtain highly accurate (<0.5%) molecular volumes VL;56 modern simulations typically match the experimental VL to within 2%. The material moduli may also provide tests of simulations; as noted in the first part of this lecture, the bending modulus KC is probably not well enough determined, but perhaps values of the area modulus KA are robust enough.19
Compared to simply running NPT simulations and comparing to the raw data, what I am suggesting in this lecture is that more insight will be obtained by running NPAT or NPγT simulations for a sequence of A or γ values. Such a sequence will allow one to estimate the area ANMR that best fits the NMR data, the area Axray that best fits the X-ray data63 and the area Aneutron that best fits the neutron data, as indicated in Fig. 5. Then, the most important validation of a simulation methodology is that Aneutron = AX-ray = ANMR, because this requires the simulation to obtain correctly three important structural features, namely, the water distribution (DB), the headgroup thickness (DHH) and the chain conformations (order parameters, NMR and also wide angle X-ray41). These are the properties that depend very sensitively on the value of A, unlike lipid volume or area compressibility KA. Of course, it would also be desirable that the common best area A would occur with zero tension, as in NPT simulations, but due to the necessity of using highly approximate water force fields, I believe that this should not be the main criterion for judging force fields or for deciding whether a simulation can be used to provide the kind of information shown in Fig. 2. Much more important is matching the different kinds of experimental data with a value of A =Aneutron = AX-ray = ANMR for each lipid bilayer system. While it is likely that the value of the effective surface tension that has to be applied will be different for different lipids, as well as for lipid bilayers with additives like cholesterol, I think it is sufficient, and may prove challenging enough, for lipid force fields to be able to achieve this limited goal for the many different lipid bilayers that have different values of A. If this can be achieved, then those best fitting values of A would be the best estimates of the area and much else besides. This idea has been suggested previously, but only X-ray data were used.65 Of course, with increasing computer power allowing better water models to be used, the ultimate goal is that NPT simulations will achieve equality for all the areas as in eqn (5).
With the exception of the effect on the bending modulus of adding cholesterol in section 2.3, so far I have focused on simple lipid bilayers. Since biomembranes consist of mixtures of many different lipids and proteins, studies of more complex model systems are of high interest. A general focus is the long standing one of lipid–protein interactions. Another general focus is membrane lateral heterogeneity, often characterized by the buzzword, rafts.
While selective deuteration in neutron scattering experiments can add more data, generally, the increase in obtainable structural data doesn’t keep up with the increase in the complexity of the model systems. Therefore, the relative value of simulations becomes potentially even greater than for simple lipid bilayers. This provides even greater motivation to obtain reliable force fields from careful studies of single component lipid bilayers. Because of the long time for lateral equilibration in atomistic simulations and the need to study systems with more lipids, coarse grained (CG) simulations have become more popular.7,49,53,54,83
I will not attempt to review the voluminous raft literature, nor will I address the controversial issues of exactly what is a raft, what are the size and time scales, whether they are primarily lipid rafts or if they require protein, or whether they are non-equilibrium structures induced by initial metabolic conditions. Whatever the answers ultimately will be to those questions, it is clear that equilibrium measurements on optically observable domains in model mixtures will be valuable to characterize molecular interactions and that those same interactions will be present in whatever biological rafts turn out to be. It is therefore certainly appropriate that this general theme of lateral organization is included in several papers in this Discussion.48–50,53,54,84,85
Fig. 7 shows a ternary phase diagram of the canonical raft mixture that resulted from a combination of measurements.86 (Note that quite different diagrams for other ternary mixtures occur in this Discussion.87) The most famous of the canonical raft diagrams was built from fluorescence imaging and NMR data.88 Recently, we added X-ray data that modified the boundaries of the three-phase and low cholesterol regions; we also made sure that our diagram conformed to Schreinemakers’ equilibrium thermodynamic rules that have often been violated by previous diagrams.86 Most importantly, we obtained the slopes of the two-phase tie lines with even greater accuracy86 than had earlier been obtained using NMR.89 Even though the construction of the full phase diagram is challenging experimentally, as well as for theory and simulations, any simulation or theory that exhibits two-phase coexistence, even for just one composition, will automatically obtain the slope of a tie line that should then be compared to experiment to help evaluate and modify the proposed interaction parameters.
Another category of more complex model systems is the addition of peptides or other membrane active molecules to lipid bilayers.3,48,52,84,85,90,91 There are two broad foci for the information one would like to obtain in such systems. One focus addresses how the additive resides in the bilayer. Is it in the hydrocarbon core or in the interfacial headgroup region or attached peripherally; what is the conformation of the additive and what is its orientation; does it form pores? Another focus asks what are the effects of the additive on the bilayer. Does the bilayer become thinner; does it become more flexible; does the additive induce local spontaneous curvature or disrupt the bilayer?
Again, I will not try to review the literature, but just mention some of the topics our lab has also worked on. Fig. 8 shows an X-ray scattering pattern before and after alamethicin was added to a bilayer. Regarding the structure of alamethicin, the additional scattering has been shown to be due to pore formation and the dimensions of the pore have been deduced consistent with of an average of about 6 monomers/pore in DOPC bilayers.92,93 Consistent with a hydrophobic matching energy penalty, the number of monomers/pore increases with thicker bilayers. Regarding the effect of alamethicin on the bilayer, there is little change in the thickness of thin bilayers, but a significant thinning of thick bilayers.75 This is also consistent with hydrophobic matching; alamethicin can tilt with no change of the thickness of thin bilayers, but it drives thick bilayers to become thinner; the experimental break point between thick and thin regimes gave 27–28 Å for the hydrophobic thickness of alamethicin. This alamethicin study75 and a recent study that added bioflavinoids to bilayers63 are examples that synergistically combine simulations and experiment.
In another study of the fusion peptide of the HIV-1 virus, we found that FP23 penetrates into the hydrocarbon core of bilayers, but does not appear to be transmembranous.64 The most striking result is that FP23 greatly decreases the bending modulus KC;94 we believe this makes it easier for the viral envelope membrane to fuse with the target T-cell membrane, forming a pore through which the T-cell becomes infected with viral RNA. This is a rather different picture than that of a three-legged space invader module punching holes using rigid spikes through the rigid skin of a target space ship; rather, ours is a soft condensed matter picture with thermally activated intermediates, which is a picture drawn from the lipid fusion literature that is more likely to appeal to physical chemists than to biochemists.
An introduction to this Discussion, and to this field, would be remiss not to emphasize the many shapes and phases obtained by lipids95 that are included in many papers in this Discussion.3,6,8,47,51,83,87,90,96 Variations from the fluid lipid bilayer phase, such as the gel phase, the inverted hexagonal II phase and the cubic phases, are relevant to the discussion of biological processes, and they are fascinating as phenomena for physical chemical modeling. Even the relatively biologically irrelevant, but physically interesting, gel phase continues to play a role in that the phase transition into the fluid phase is an easy measure of the effect of additives90 and this phase is a valuable testing ground for techniques.8,47 The hexagonal II phase83 and the cubic phase6 are used in discussing membrane fusion.97,98 Other structures, such as tubules3,51 or fibers,96 are relevant biologically in that membranes undergo morphological transformation required for cell function. This great diversity of phases and shapes provides a huge database for further physical chemistry study.
I thank Norbert Kučerka for maintaining the SIMtoEXP website and Richard Pastor, Rick Venable, Jeff Klauda, Scott Feller and Jonathan Sachs for transmitting data from published simulations in an easy to use SIMtoEXP format. I also thank Olle Edholm and especially Anthony Braun in the Sachs’ group for running specific, as yet unpublished, simulations at my request. Errors of interpretation herein of any of these published or unpublished simulations are, however, my responsibility. Research reported in this publication was supported by the National Institute of General Medical Sciences of the National Institutes of Health under award number R01GM044976. The content is solely my responsibility and does not necessarily represent the official views of the National Institutes of Health.
†It is interesting that my method for interpreting the NMR data,81 which ignored reversals in direction in the hydrocarbon chains subsequently found both experimentally82 and in simulations (indicated in the terminal methyl distribution in Fig. 2), nevertheless obtained the value A = 62 ± 2 Å2 for DPPC at 50 C, which agrees well with the most recent value58,59 of A = 63 Å2. However, the NMR data are less valuable when there are unsaturated chains, and there are additional issues74 that lead me to conclude that NMR should also not be the primary way to estimate A.