Method for obtaining structure and interactions from oriented lipid bilayers

I. INTRODUCTION

Lipid bilayers have been much studied as the prototypical biomembrane. One focus of such studies is the structure of lipid bilayers and how that structure varies with the thermo-dynamic phase and with the chemistry of different lipids [1,2]. Another focus is on the interactions between two bilayers [3,2] and on the mechanical properties of bilayers [4].

For these purposes many groups have used x-ray scattering from lipid bilayers in the model system of multilamellar vesicles (MLVs) [1-3,5]. Each MLV consists of bilayers arranged (topologically at least) as more than 200 (usually many more) concentric spheres with an average repeat spacing D of order 60 Å , of which about 25 Å is water, depending upon the lipid and the phase. For scattering purposes MLV samples effectively consist of independent domains, each of which is a stack of parallel bilayers, but the normals to the bilayers in the independent stacks are isotropically distributed in space. Such powder samples diffract relatively weakly and the intensity falls off rapidly with q because of the Lorentz factor which scales as 1/q². For samples in the fully hydrated, most biologically relevant, fluid (L_α) phase, often only h=2 orders of diffraction are observable (q_z,_max ≈0.2 Å⁻¹). Dehydrating the sample by applying osmotic pressure increases the number of observable orders, but there is an osmotic pressure correction to structure that must be applied [2,6,7]. Furthermore, for at least one lipid (DOPC), the structure changes drastically upon mild dehydration (relative humidity near 96% [8]; for lower relative humidity, diffraction data are not relevant for the desired bilayer structure [5].

The intensity of diffraction is much greater from stacks of bilayers that are oriented with their normals aligned along one axis. However, until quite recently there has been concern about studying oriented samples by x-ray diffraction because they could not be hydrated from water vapor to the full D spacing of samples immersed in water. This “vapor pressure paradox” has recently been shown to be an experimental artifact [9,10] and sample chambers suitable for x-ray studies have now been built that obtain full hydration [11]. Figure 1 shows the kind of data that hold the promise of providing far more information than the MLV samples. The purpose of this paper is to show how data like those in Fig. 1 can be analyzed.

The first reason that the data in Fig. 1 are so promising is that there is intensity for a q_z range at least up to 0.6 Å ⁻¹, more than twice as far as for the comparable MLV samples [12]. The second reason is that the deep minima in the intensities near q_z=0.28 and 0.44 Å⁻¹ immediately show where there are zeros in the single bilayer form factor F(q_z), defined by

F(q_z) provides information about the structure of a bilayer centered at z=0 with average normal direction z through the contrast in the electron density ρ_I*(z) relative to that of water ρ_W*. Third, in addition to data in the q_z direction, there is clearly distinguishable information in Fig. 1 along the q_r direction which, for MLV samples, is lost by convolution into a powder-averaged q dependence [13]. Finally, the data are not confined just to the vicinity of the peaks as in our previous studies on MLVs [13], although other groups [14-16] have obtained and analyzed intensity data that are continuous in q for MLV samples.

This work is closest in many ways to the study of Lei et al. [17] and has benefited from it (see especially Lei [18)]. There are, however, substantial differences which require a separate development. Lei et al. studied a surfactant system where the bending modulus could be made small by adding cosurfactant so that their inverted bilayers could swell to large D values. This took the system into the hard confinement regime where the Helfrich theory of interactions [19] is appropriate, whereas our systems have much smaller D and are in the soft confinement regime [20,21]. Because fluctuations were large, at most two peaks were observed and data extended only to about 0.25 Å⁻¹. Structure was not a goal, so the electron density of the bilayer ρ_I*(z) was assumed to be a constant which is an adequate assumption for their q_z range, but not for ours for which the extended q range provides more information about the structurally useful modulation of ρ_I*(z) across the bilayer. Their major goal, which is also one of ours, was to obtain the bending modulus K_c and the bulk (interaction) modulus B. These two moduli are related to two parameters that were obtained from the scattering data. The first of these is the Caillé parameter [22]

The second, which involves B and K_c in a different combination, is an in-plane correlation length [23,17]

which is related to the de Gennes penetration length λ[24] by ξ²=λD. The Caillé parameter has been obtained in many studies, including MLVs, by analyzing the tails of the diffraction peaks [13,25]. Obtaining ξ is even more challenging [26]; in particular, powder data are very insensitive to ξ[13]. Lei et al. [17] obtained ξ by analyzing the intensity in the very small angle regime (q_z<0.03 Å⁻¹). However, as they noted, this small-angle scattering becomes weak in the more weakly fluctuating soft confinement regime appropriate to our lipid bilayers. We therefore did not attempt to obtain data at such small angles (this small-angle regime is obscured by the beam block in the data in Fig. 1). Instead, we show how both η₁ and ξ (equivalently, B and K_c) can be obtained from the q_r dependence of the data at large values of q_z .

Although the data in Fig. 1 contain a great deal of information, including the structural information in F(q_z) and the interaction information in η₁ and ξ, the data are also affected by the domain size. Following [13], we define L_z to be the average domain size in the direction of the average bilayer normals. For cylindrical domain shapes, L_r is defined to be the average radius in the transverse direction along the plane of the bilayers. We also consider distributions in these two sizes. The data also depend upon inevitable experimental artifacts such as x-ray beam size, mosaic spread, and instrumental resolution. We will initially ignore the experimental issues, except to note that the instrumental resolution for the data in Fig. 1 was limited by the pixel size of the charge-coupled device (CCD) detector, about 0.001 Å⁻¹, which is 10 times coarser than our best resolution for MLV samples [12]. Therefore, our analysis must not rely upon the most intense data very near sharp diffraction peaks, but must instead focus upon the weaker intensity outside the peak regions.

Since the data are affected by so many parameters, the analysis should ultimately consist of a simultaneous fit of a model to the global data. However, in this paper we prefer to establish the feasibility that such a program even has a chance to succeed by identifying different regions of the data from which a smaller subset of parameters can be estimated separately. That this piecemeal approach might work is indicated by the fact that the high-q_z data are so strongly affected by fluctuations that all the scattering is diffuse whereas the smaller-q_z data exhibit sharper and narrower peaks, more typical of ordered samples. Another indication is that the intensity I(q) is a product

of a structure factor S(q), which contains the B and K_c information, and a form factor squared, which does not have any q_r dependence, as well as the customary q_z⁻¹ Lorentz factor for oriented samples. This implies that the q_r dependence must therefore come from S(q). Therefore, the determination of the structural parameters in F can be separated from the remaining parameters involved in S.

In Sec. II the basic model of discrete bilayer stacks is reviewed and it is shown how the correlation functions that affect x-ray scattering are calculated. It is shown how the structure factor S(q) is calculated in Sec. III. Regions in which the domain sizes and distribution have little effect on S(q) are identified in Sec. IV. It is then shown in Sec. V that both η₁ and ξ can be determined from data in these regions. We return in Sec. VI to the other regions where the domain sizes do affect S(q) and show how the sizes can be estimated from data in these regions. It is then shown in Sec. VII how the previous sections enable determination of the form factor F(q_z). Finally, in Sec. VIII this new analysis is used to provide results for the data in Fig. 1.

II. CORRELATION FUNCTIONS FOR FLUCTUATIONS

The system consists of n=1,...,N bilayers stacked in the z direction. The center of each bilayer has an average position z_n=nD and each bilayer extends over a diameter L_r in the r=(x,y) plane. Fluctuations in the stacking of the bilayers are described by the quantities u_n(r) which are spatial deviations in the z direction of the center of the nth bilayer from its average position at location r. The original Caillé theory [22] employed a continuum model of smectic liquid crystals which has gradually been replaced by the following, more realistic, discrete free energy functional [27,20,17,28,21]

where the first term accounts for the bending energy of individual bilayers and the second term is a harmonic approximation to the interactions between adjacent bilayers [20,28], parametrized by the modulus B which has units of energy/length ⁴. Our B is related to the Caillé three-dimensional bulk modulus B₃ by B=B₃ /D.

To treat x-ray scattering it is necessary to calculate the height-height pair correlation function

which has no azimuthal dependence because the bilayers are in-plane fluids. We will avoid the ambiguity that arises from imposing boundary conditions on a finite-size domain by calculating the correlation functions for an infinite domain and assume translational invariance. The restriction to finite-size domains will be made in Sec. III by summing the infinite domain correlation functions only over pairs that are contained in such a domain. Then, for the discrete model in Eq. (5) it has been shown [17,28] that

where J₀ is the zeroth-order Bessel function. Defining x =(ξk_r)²/2 and integrating over k_z using the calculus of residues [18] yields

where q₁=2π/D is the position of the first diffraction peak. Equation (8) shows that η₁ is just a scaling factor for δu_n(r). This is computationally convenient because a table of correlation functions can be calculated once to use with different values of η₁. Something similar occurs with the parameter ξ , but this is somewhat less obvious, since ξ appears both inside the integral and in the upper limit in Eq. (8). Fortunately, the upper limit of integration can be extended to ∞ with high accuracy [18], as we have verified. To remove ξ from inside the integral, we change variables to t =r/ξ so that the height-height correlation functions finally can be written

in which ξ appears only as a scaling factor of r. This means that we can calculate a table for δu_n(r) only once for certain values of r=ξt, n, η₁, and ξ . Then, we can use this table for other values of η₁ by simply rescaling the values of δu_n(r) by a constant factor. And we can use this table of δu_n(r) for other values of ξ by rescaling the r dimension. We calculate δu_n(r) in Eq. (9) numerically as a function of r and n on a grid and save the values in a table. (It is computationally useful to use logarithmic steps in r.) We then perform smooth polynomial interpolations for δu_n(r) for any r between the grid points.

Caillé [22] proposed an approximate correlation function for r>ξ,

where γ is Euler’s constant and E₁(t) is the exponential integral. In Fig. 2 we compare this approximation to our accurate results for n=0. Even for this worst case good agreement is obtained for r>1000 Å. For n>0 the approximation works well for even smaller values of r. However, precise numerical calculation is required in the vicinity of r=0 and n=0 and also high accuracy is needed for detailed analysis of the structure factor, so we use Eq. (9) to calculate δu_n(r) for r<2000 Å and n<30. We therefore use the Caillé approximation, Eq. (10), for r>2000 Å or n>30, having checked that the fractional difference is less than 10⁻⁴.

III. CALCULATION OF THE STRUCTURE FACTOR

The electron density ρ_n* of the nth bilayer will be written

where ρ_I*(z) is the electron density profile of the bilayer centered at z=0. As has been discussed before [13,29], there is little reason to suppose that the local fluctuations in the bilayer shape δρ_I* are correlated with the long-wavelength fluctuations u_n(r), at least when the bilayers are fully hydrated [5]. However, there is a geometric correction to the projected thickness of undulating bilayers along the average bilayer normal which is correlated with |▽u_n(r)| and, when under osmotic pressure, the local bilayer thickness is anticor-related with |▽u_n(r)| [5]. Both effects are limited to about 2% for typical lipid bilayers and will be ignored. Then, the basic scattering formula

can be written as Eq. (4) [13,17] where F(q_z) is given by Eq. (1) and the structure factor S_L(q_z ,q_r) is

where G is the scattering pair correlation function:

Since the lipid bilayer is organized as a fluid within the plane of each bilayer and there are no local correlations between molecules in adjacent bilayers, each bilayer is an in-plane powder, the azimuthal angle θ can be integrated out, and the in-plane vector r can be replaced by the in-plane distance r. As already discussed in Sec. II, the correlation functions depend only upon (n—m) and (r—r’), so the structure factor for a domain of size L_r and L_z=ND becomes

where F_r(x), which is a finite-size effect function in the r direction [18], is zero if x>1 and for x≤1

We next assume that the distributions of domain sizes L_r and L_z are independent, so that

where P_z(L_z) and P_r(L_r) are distribution functions for L_z and L_r . (One might alternatively assume that the L_r and L_z sizes are correlated in some fashion.) We will use the Gaussian distributions

Equation (17) can be simplified by introducing effective finite-size factors

Thus, the final equation for S(q_z ,q_r) becomes

To calculate S(q_z ,q_r) from Eq. (22) numerically we first calculate G using the table for δu_n(r) in Sec. II. We also calculate H_z(z) and H_r(r) using Eqs. (20) and (21) and save them in tables. We then calculate the integral and sum in Eq. (22) numerically. The accuracy of our calculations has been checked in two ways. The same integration was performed using different tables calculated using a finer grid for r in δu_n(r) and different integration procedures of the Romberg type have been used, with cutoffs adjusted to provide an accuracy of 10⁻⁴.

Figure 3 shows calculated structure factors where the values of η₁ and ξ were chosen to correspond to K_c=10⁻¹² ergs and B=2.8×10¹² ergs/cm⁴, which are typical for lipid bilayers [28]. As expected [17,13], the intensities of the higher orders rapidly decrease and almost disappear by h=5 even at q_r=0. As q_r increases, the higher-order peaks disappear completely, but the first two peaks are seen even at q_r=0.04 Å⁻¹. Curiously, the shape of the peaks as a function of q_z becomes narrower and more Bragg like as q_r increases. In the next several sections we will discuss how changes in the various parameters affect the structure factor shown in Fig. 3.

IV. EFFECT OF DOMAIN SIZES

Figure 4 shows how the structure factor is affected by the average domain size in the z direction which henceforth will be designated just as L_z with no overbar. For all comparisons we normalize S(q_z ,q_r) by the amount of material in the bilayer stack. This normalization gives the same values for S(q_z ,q_r) for large values of q_z . As expected, the peaks also become higher with larger domain size. Because of the overall conservation of scattering intensity, the diffuse scattering between peaks becomes lower with an increase in L_z . Most importantly, the dependence on L_z becomes smaller as q_z and/or q_r increases. There is almost no dependence between the h=3 and h=4 peaks when q_r=0. For q_r>0.01 Å⁻¹ there is almost no dependence on L_z except very close to q_z=2πh/D for h=1 and h=2.

Figure 5 shows that L_r affects the structure factor quite differently from how L_z affects it. An increase in L_r results in an increase of the whole S(q_z) curve when q_r=0, more so at low q_z values than at high q_z . This does not violate conservation of total scattering intensity since an increase in L_r makes the peaks in S(q_r) sharper and higher, therefore reducing the intensity in the q_r direction. Most importantly, and similar to the effect of L_z above, the dependence on L_r becomes smaller as q_z and/or q_r increase. There is almost no dependence between the h=3 and h=4 peaks when q_r=0. For q_r>0.01 Å⁻¹ there is almost no dependence on L_r except very close to q_z=2πh/D for h=1 and h=2.

V. DETERMINATION OF η₁ AND ξ

The results in the previous section show that the domain size parameters can be ignored in large regions of the data where the fluctuations and diffuse scattering dominate. The important unknowns in these regions are therefore η₁ , ξ , and the form factor F(q_z). Since F depends only upon q_z [Eq. (4)], we now concentrate on the q_r dependence of the structure factor to evaluate the possibility of finding η₁ and ξ without having to do a global fit that includes F(q_z) and the domain size parameters.

Figure 6 shows S as a function of q_r in regions where the dependence on domain sizes is negligible. Since the overall amplitude of these curves would be set by the as yet unknown form factor, the curves are normalized at q_r=0.01 Å⁻¹ in order to concentrate on the q_r functional dependence for larger q_r . Figure 6 shows that both parameters η₁ and ξ affect S(q_r). There is an important difference, however, between the effect in Fig. 6(a), which is for a q_z slice (constant q_z and varying q_r) at the h=fourth-order peak, compared to the effect in Fig. 6(b), which is for a q_z slice that goes midway between the third and fourth peaks. In Fig. 6(a) the effect of changing η₁ by a factor of 2 is relatively larger than changing ξ by the same factor and the reverse is true in Fig. 6(b). Therefore, the effect of changing the two parameters is not functionally degenerate and it is therefore plausible that both parameters can be extracted from two q_z slices of data.

Programs were written to fit data from two or more q_z slices and the programs were tested on calculated data with random noise comparable to the noise in the data in Fig. 1. (We note that data can be taken in future with less noise.) Two slices of calculated data are shown in Fig. 7. The fitted curves are also shown in Fig. 7 and the fitted values of η₁ and ξ agree well with the input values.

We note that a choice of two q_z slices should include at least one with q_z near a peak value hq₁ and one from between peaks, such as (h−½)q₁, as shown in Figs. 6 and 7. Choosing two slices, both between successive peaks, does not give sufficient contrast to obtain good parameter values when noise is included. We have used more than two q_z slices to increase the accuracy of the determination of the parameters, but two slices suffice to illustrate the principle.

VI. DETERMINATION OF DOMAIN SIZES

In Sec. IV the emphasis on the effect of domain sizes was to find the scattering regions where the domain sizes L_z and L_r had little effect so that the evaluation of η₁ and ξ could be carried out independently as shown in Sec. V. We now turn to the scattering regions near and between the low-order peaks where the domain size has a large effect. This will give us a method to estimate domain sizes. This procedure uses the values of η₁ and ξ that will now be assumed to be determined by the method of Sec. V. We cannot assume the values of the form factors F(q_z), so the method in this section will again focus on the q_r dependence of the structure factor but now at low values of q_z .

Figure 8(a) shows the dependence of S(q_r) on L_z for two q_z slices, one that goes through the h=1 peak and one that goes halfway between the h=1 and h=2 peaks. We have chosen to normalize these curves to unity at q_r=0.01 Å⁻¹, which is the same value of q_r used for normalization in Sec. V, since there is little dependence of S on domain size for q_r>0.01. With this normalization the dependence of S on L_z near q_r=0 is opposite for the two q_z slices. Figure 8(b) shows a much different dependence of S(q_r) on L_r for the same normalization. For both q_z slices an increase in L_r results in an increase of S near q_r=0.

The aforementioned difference in the effect of domain sizes L_z versus L_r on the behavior of S(q_r=0)/S(q_r =0.01) at different q_z values indicates that these two parameters are functionally inequivalent and therefore both L_r and L_z should be extractable from the data, as we have verified for calculated data in which we have varied L_z and L_r . However, this assumes that we know the functional form of domain sizes. Our experience with powder samples has led us from Gaussian distributions [12] to exponential distributions [28]. It is also not clear how L_r and L_z should be correlated. We think it is unlikely that these kinds of issues will be definitively resolved. Fortunately, these domain details are also the least important from the fundamental point of view, since, due to the exigencies of sample preparation, they are strongly variable even for nominally identical samples, as was earlier emphasized for powder samples [25]. (One would expect oriented samples to have even more variation.) Therefore, we are content to obtain a rough estimate of domain sizes and distributions and to focus on data in the regions where details of the domains are not important.

VII. DETERMINATION OF THE FORM FACTOR

Since S(q) is determined from the methods in the previous sections, we will now use Eq. (4) to calculate the form factor F(q_z), which is the only remaining unknown. As noted at the end of the previous section, because details of domain sizes and distributions are the most problematical quantities to obtain precisely, it is appropriate to determine F(q_z) by analyzing data in regions in q where such details are unimportant. We have chosen to use I(q) data for the region 0.035<q_r<0.040 Å⁻¹ where intensity data are also strong enough to perform reasonable fits. For each value of q_z , F(q_z) is chosen to provide the best overall fit to Eq. (4) in this q_r range.

We have tested this method of finding F(q_z) using calculated data with added noise. As expected, the method works well when one knows η₁ , ξ , L_z, and L_r, precisely. Figure 9 shows the result of a more challenging test which supposes that an incorrect value of L_r and L_z was obtained from the method in Sec. VI. In this test the correct value of the form factor is F=1 for all q_z . Figure 9(a) shows that the estimated F(q_z) is quite good everywhere except very near the diffraction peaks, even though L_z was misset by a factor of 5. Figure 9(b) indicates that missetting L_r has a larger, more systematic effect, but it is still confined to the vicinity of each diffraction peak. We note that the values of F in Fig. 9(b) are consistently a few percent too high for most values of q_z not near a diffraction peak, but such a uniform error is unimportant since experimental intensities are relative. Since F(q_z) is a smooth function on the q_z scale of each of the deviations in Fig. 9, one can easily interpolate accurate values of F(q_z) for these regions where one expects the method to be inaccurate.

VIII. APPLICATION TO DOPC DATA

We have applied the analysis methods developed above to the data shown in Fig. 1. The instrumental resolution of these data was limited to δq=0.001 Å⁻¹ by the pixel size (50 μm of the CCD detector and its placement (22 cm) from the sample. The sample volume was effectively 10 μm ×500 μm and the latter length was convolved into the calculations to fit the data. The small mosaic spread, estimated to be less than one degree, was not included in the fit.

The results for the fluctuation parameters are η₁=0.08 ±0.01 and ξ=59±2 Å. The value of η₁ may be compared with the value η₁=0.10±0.01 obtained from the analysis of peak shapes of unoriented fully hydrated MLV samples [6]. Using Eqs. (2) and (3) gives K_c=0.73±0.04×10⁻¹² erg and B=6.0±0.7×10¹² ergs/cm⁴. This value of K_c may be compared with a recent value K_c=0.85±0.1×10⁻¹² erg obtained using the aspiration pipette method on giant unilamellar vesicles at the lower temperature of 18 °C [4]. The value of B has not been previously obtained for fully hydrated DOPC, but the new value generally agrees with the range of values reported for egg lecithin [28].

The results for domain sizes are L_z=0.18 μm and L_r=0.5 μm. The thickness of the sample on the curved glass substrate is at most 10 μm, so the value of L_z appears reasonable. We expected a larger size for L_r . It is noteworthy, however, that L_r is still considerably larger than ξ and L_z is considerably larger than the interbilayer water spacing of 0.0025 μm, so a theory of long-range fluctuations is appropriate for this sample.

The result for the form factor is shown in Fig. 10. It is compared to a published form factor obtained for partially dehydrated MLV samples which was put on an absolute scale [6]. The overall scale factor for the new form factor from the oriented sample is set to match the previous form factor in the first lobe region q_z<0.3 Å⁻¹ where the MLV results are most accurate. The new results become very noisy for q_z<0.15 Å⁻¹ due to uncertainties in subtracting background which is higher for small q_z due to air scattering from the beam. The older results in the region of the third lobe (q_z>0.44 Å⁻¹) were obtained with very few data with large statistical deviations [6]. The new results for F appear to be more accurate in the high-q_z region.

IX. CONCLUSIONS

It has always been attractive to study lipid bilayers stacked in oriented samples because they scatter more strongly. Now that the vapor pressure paradox has been over-come experimentally [9,10], there is added motivation because these samples can be fully hydrated in the biologically relevant fluid phase. However, the intensity is affected by a convolution of many factors. The analysis presented here shows how to extract both the parameters K_c and B from x-ray scattering data from such samples in the fluctuating L_α phase. These parameters appear in the fluctuational free energy that has previously been shown [12,28,6] to describe powder samples; however, only one of the parameters could be obtained from the powder data. Being able to obtain both parameters will enhance future studies of the basic interactions between lipid bilayers in which the osmotic pressure will be varied systematically to obtain B as a function of the distance between bilayers [3,2,28].

Our analysis shows how to obtain the continuous form factor F(q_z) in a way that is not very sensitive to uncertainties in determining domain size. Previous determinations of the continuous form factor using powder samples relied on only a few orders from many samples, each with a different D. In contrast, the results in Fig. 10 come from just one sample. The relative wealth of q-space data in Fig. 10 raises several possibilities for obtaining electron density profiles for lipid bilayers. In a subsequent paper we will consider some of these possibilities, employing F(q_z) results that we hope will extend to even higher q_z .

It has also not escaped our attention that this method should be applicable to more complex biomembranes that include mixtures of lipids, sterols, or antibiotic peptides in the biologically relevant, fully hydrated state.