INTRODUCTION

TOP ABSTRACT INTRODUCTION METHYLENE TRAVEL AVERAGE CHAIN LENGTH HYDROCARBON THICKNESS AREA PER MOLECULE DISCUSSION REFERENCES

Simple lipid bilayers have been a challenge for quantitative structure determination. For example, experimental values for the average area per lipid A for fully hydrated dipalmitoylphosphatidylcholine (DPPC) at 50°C have ranged from 56 to 72 Å² (Nagle, 1993), with corresponding uncertainty in the average thickness of the hydrocarbon core, 2D_C. The same range of uncertainty has come from the two primary experimental methods, namely x-ray diffraction and NMR. Recently, primary emphasis has been on the x-ray technique, where A_DPPC = 62.9 ± 1.3 Å² was obtained (Nagle et al., 1996). This result was in agreement with a reanalysis of NMR order parameter data, which gave A_DPPC = 62 ± 2 Å² (Nagle, 1993). However, it was emphasized in that reanalysis that there were assumptions that could and should be tested with simulations. This paper carries forward that program.

Simulations give a much more detailed view than any experiments on lipid bilayers (Tobias et al., 1997; Tieleman et al., 1997). However, given the uncertainties in force fields and the restriction to nanosecond time scales, one should not necessarily expect simulations to obtain accurate values for all quantities of interest. The use of simulations that we envision is to test relations between simulated microscopic quantities that cannot be obtained from experiment. Such relations are commonly used to obtain quantities of interest from raw data. Simulations can also inspire new relations, as we show in this paper. This use of simulations does not require that all of the force fields be exactly correct or that the simulated quantities of experimental interest agree perfectly with experiment. The interplay of simulations that is envisioned is that simulations will help guide experimental analysis, which will then provide more accurate quantities of interest, which will then help tune the force fields used in the simulations.

The main relation studied in this paper is between the average order parameters S_n and the average travel per methylene D_n along the normal to the bilayer. The conventional relation (Schindler and Seelig, 1975; Thurmond et al., 1991; Nagle, 1993) can essentially be written as

(1)

where D_M is the maximum possible travel. A different relation has also been used (De Young and Dill, 1988):

(2)

However, a recent simulation (Berger et al., 1997) shows (their figure 7) that neither relation is especially good, and this paper reports essentially the same result for a different simulation. Even though the values of A that were obtained in their simulation using the method of Nagle (1993) were in fairly good agreement with the actual simulated A, it is important to elucidate the errors, including fortuitous compensations.

The simulation we analyze in this paper has been described previously (Tu et al., 1995). Briefly, the simulation was performed at constant number, pressure, and temperature (NPT), with N = 32 DPPC lipids per monolayer and n_W = 28 waters/lipid, P = 0, and T = 50°C, using a Nose-Hoover thermostat. The SPC/E water potential and an all-atom description of the lipids were used in the CHARMM program, with Ewald sums for long-range interactions. Since publication (Tu et al., 1995) the simulation has been continued to 2 ns.

	METHYLENE TRAVEL

TOP ABSTRACT INTRODUCTION METHYLENE TRAVEL AVERAGE CHAIN LENGTH HYDROCARBON THICKNESS AREA PER MOLECULE DISCUSSION REFERENCES

The first main result is illustrated in Fig. 1, which plots the conventional CD₂ order parameters S_n versus the travel of the nth methylene along the bilayer normal. Each simulated data point in this figure gives molecular dynamical averages for one specific methylene, with the different points corresponding to different carbon numbers n, different chains (sn-1 and sn-2), and different monolayers (upper and lower) in the simulated bilayer. The solid straight line shows that the conventional prediction (Nagle, 1993) works fairly well for values of the order parameter S_n  0.2 that are close to the experimental plateau region, but the slope of the data is clearly smaller than the conventional slope. This point has previously been made for a different simulation by Berger et al. (1997); they presented essentially this type of figure (except that they plotted the molecular order parameter instead of S_n) with very similar results. Our first main result is the derivation of a new formula that gives the result shown by the solid curve in Fig. 1.

View larger version (20K): [in this window] [in a new window]   FIGURE 1   Symbols: simulated average order parameters Sn versus average molecular travel Dn (in Å). Curved solid line: New Eq. 24. Straight solid line: conventional Eq. 1. Dashed line: Eq. 2. Dotted line: Eq. 14.

The derivation of the new formula was inspired by detailed results of the simulation that have not previously been presented. As with any fluctuating statistical mechanical ensemble, one should consider distribution functions. A fairly general one that we have analyzed for the hydrocarbon chains is denoted p_n(z, D). The discrete index n is the carbon number, which goes from 2 for the methylene next to the carbonyl to 16 for the terminal methyl. The continuous variable z is the distance along the bilayer normal, where z = 0 is the center of the bilayer and with the positive sign for the direction toward the headgroup region for that chain. The instantaneous position of carbon n is denoted by z_n. The continuous variable D is the methylenic travel; for carbon number n the instantaneous travel D_n is defined by

(3)

as illustrated in Fig. 2. From the distribution function p_n(z, D) the average position of carbon n is given as

(4)

Fig. 3 shows that z_n has upward curvature. This corresponds to D_n decreasing with n, which conforms to the expectation that the chains become more disordered near the terminal methyl end.

View larger version (22K): [in this window] [in a new window]   FIGURE 2   For the nth CH2 group, D defines the travel along the bilayer normal z,  defines the local methylenic tilt, and  defines the rotation about the local methylene axis z'. Note that the local axis is generally different for each carbon n.

View larger version (12K): [in this window] [in a new window]   FIGURE 3   Average carbon position zn (solid squares, in Å) for all chains.

For many purposes it suffices to consider only the reduced distribution function

(5)

For example, the average travel of carbon n, irrespective of where it is, is then given by

(6)

It is also important to consider the mean square travel,

(7)

These quantities are shown in Fig. 4.

View larger version (13K): [in this window] [in a new window]   FIGURE 4   Simulated data for n () and Dn () averaged for all chains.

To consider averages of the order parameter S_n, it is necessary to consider another variable in the probability distribution functions. Fig. 2 shows the local axis z' that goes through carbons n  1 and n + 1 and makes the angle  with the bilayer normal along z. Rotation around the local axis z' is measured by the angle . A straightforward calculation (see Eq. 3a; Nagle, 1993) gives the combined order parameter S_n of both deuterons on the nth methylene corresponding to this orientation:

(8)

where cos  = D_n/D_M involves only the previous variable D_n and the maximum travel along the z axis, which is D_M = 2.54 Å for undistorted saturated hydrocarbon chains. The additional variable is sin², for which average values for narrow ranges of D are given by

(9)

Fig. 5 shows values of sin²_n(D) for two values of n for the upper and lower monolayers separately. The fact that sin²_n(D) is closer to 0.4 for the upper monolayer and closer to 0.5 for the lower monolayer indicates incomplete thermal equilibration. This variable was considered by van der Ploeg and Berendsen (1983), with the suggestion that sin²_n(D) was less than 0.5. 

View larger version (17K): [in this window] [in a new window]   FIGURE 5   Simulated sin2n(D) as a function of D for n = 6 (, ) and n = 12 (, ) averaged over sn-1 and sn-2 chains. , , The upper monolayer; , , the lower monolayer.

After integration of Eq. 8 over , the average order parameter

(10)

becomes

(11)

Fig. 6 shows simulation results that have been binned according to the values of D, for both the upper monolayer and the lower monolayer. Fig. 6 also shows theoretical curves for S_n(D) that were obtained from Eq. 11 for different values of sin²_n(D). Clearly, the curve using a random distribution sin²_n(D) = 0.5 fits quite well for 1.5 Å < D < D_M, where the counting statistics are good (see p(D) curves in Fig. 6), and even reasonably well for smaller values of D, where the statistics are much poorer, because there are few groups with large negative values of D; such groups have been called upturns (Nagle, 1993). (The occurrence of values of D greater than D_M is due to molecular distortions caused by thermal fluctuations and intermolecular interactions.)

View larger version (20K): [in this window] [in a new window]   FIGURE 6   Simulated S(D) averaged over all carbons in the lower monolayer () and the upper monolayer () are compared to Eq. 11. , the result for sin2(D) = 0.50; - - -, the results for 0.46 and 0.54. The curves labeled p(D) show the probability distribution function for the lower monolayer () and the upper monolayer (- - -).

Based on Fig. 6, we now adopt the approximation sin²_n(D) = 0.5. Then,

(12)

Further integration over D yields the usual average order parameters,

(13)

The consistency of Eq. 13 is tested in Fig. 7, which shows that the results for the direct simulation of S_n are very close to the values of S_n obtained from Eq. 13 and direct simulation of _n.

View larger version (13K): [in this window] [in a new window]   FIGURE 7   Comparison of Sn obtained from Eq. 13 (+) with directly simulated Sn ().

Equation 12 emphasizes that S_n(D) is a quadratic function of D, so that S_n is a function of D_n² = _n² rather than of D_n². This just reflects the well-known fundamental fact that the order parameter is a second-order Legendre polynomial, whereas the travel is first order. One might nevertheless consider approximating D_n by _n:

(14)

where the equality comes from inverting Eq. 13. This is a poor approximation, as shown first by Fig. 4, which compares D_n and _n, and, more importantly, by the lower curved line in Fig. 1. The reason for this is that _n is generally greater than D_n because p_n(D) has a nonzero width _n:

(15)

Because _n² > 0, the methylenic travel given by Eq. 14 is overestimated. This is a basic mathematical problem. Upturns, defined as groups with D < 0, are involved, in so far as they broaden the distribution and increase _n for larger n, as indicated by the increasing difference between _n and D_n in Fig. 4, but they are not the fundamental problem.

To calculate D_n from _n, we need a way to estimate the mean square deviation _n. Such a way is suggested by considering the distribution function p_n(D) in Fig. 8, which shows that, over the most significant range of D, p_n(D) behaves roughly exponentially:

(16)

The relevant parameter is the decay length _n, and the parameter _n is just a normalization factor. It may be noted that the decay length _n increases with the carbon number n in Fig. 8. Furthermore, secondary peaks in p_n(D) occur near maximally negative values of D; these are clearly deviations from the functional form in Eq. 16. These are due to upturns, which are infrequent for small n, but become more numerous, although they are still less than 10%, for large n.

View larger version (27K): [in this window] [in a new window]   FIGURE 8   Number distributions pn(D) for n = 3 (smallest decay length) to n = 15 (largest decay length).

Assuming Eq. 16, we have

(17)

(18)

With no further approximation, Eq. 17 can be solved numerically to give _n. However, analytical expressions are more appealing, and for this reason we extend the above integrals from D_M to (valid for _n D_M). Then,

(19)

(20)

(21)

Solving for _n in terms of _n yields

(22)

(23)

The square roots above are imaginary for _n² < D_M²/2, and this sets a limit to the applicability of Eq. 23. For such small values of _n² most of our previous approximations are questionable. In this limit, _n would be on the order of D_M/2, meaning a broad p_n(D), which is poorly modeled by Eq. 16. For carbons at the end of the chain we may expect more complicated probability distributions. As shown by the simulations, the maximum in p_n(D) deviates from D_M, and peaks at negative values of D start to develop because of upturns. However, for another comparison note that a constant (isotropic) distribution between D_M and D_M gives _n² = D_M²/3; this is only a little less than the applicability limit, so the restriction _n² > D_M²/2 is not so severe.

Assuming that the approximations hold, we can combine Eq. 23 with Eq. 13 to give our final expression of D_n as a function of S_n:

(24)

where the corresponding limiting point of applicability is S_n = 1/8 = 0.125. In these simulations, only S₁₅ is close to this limit, as can be seen in Fig. 1, which plots Eq. 24 and compares to direct simulations. The usual way to compare is shown in Fig. 9, which compares the direct simulation results for D_n with those obtained from Eq. 24, using the simulation results for S_n. These two results deviate for large n because the exponential form used in Eq. 16 is not accurate for methylenes near the terminal methyl, as shown in Fig. 8. Moreover, the approximation _n D_M used to obtain Eq. 19 breaks down for n > 10. The standard Eq. 1, shown in Fig. 9, deviates at both high n and low n.

View larger version (14K): [in this window] [in a new window]   FIGURE 9   Solid squares show directly simulated Dn (in Å). Estimates using Eq. 24 are shown with open circles, and the + symbols show estimates from Eq. 1.

	AVERAGE CHAIN LENGTH

TOP ABSTRACT INTRODUCTION METHYLENE TRAVEL AVERAGE CHAIN LENGTH HYDROCARBON THICKNESS AREA PER MOLECULE DISCUSSION REFERENCES

In this section we extend the results for methylenic travel to longer segments of the hydrocarbon chain. Of particular interest is the average length of the entire hydrocarbon chain. This length is a little awkward conceptually because it should include a poorly defined piece beyond the terminal methyl n = 16. A more precisely defined length will be called L^*_C, which is defined to be the average distance along the bilayer normal between the first methylene carbon and the terminal methyl carbon:

(25)

For the present simulation z₂ = 13.78 Å and z₁₆ = 1.32 Å, giving L^*_C = 12.46 Å. To estimate the average full length of the chain L_C, we first add 0.547 Å to L^*_C to account for half of the projected distance from the first methylene to the carbonyl carbon. We then add 1.5z₁₅  z₁₆ to account for the extra length of the terminal methyl; the rationale for this is that the terminal methyl volume V_CH3 is about twice as large as methylene volume V_CH2 (Nagle and Wiener, 1988; Petrache et al., 1997). The estimated chain length obtained directly from the simulations is then L_C = 13.99 Å.

To obtain chain lengths from NMR requires using either our new Eq. 24 or the conventional Eq. 1. We first report results for L^*_C, which are simply obtained by summing D_n over all odd values of n from 3 to 15. The result from Eq. 24 is L^*_C = 12.49 Å, and the result from Eq. 1 is L^*_C = 12.42 Å. Despite the deviations in Fig. 9, both approximate results from S_n agree well with the direct result, L^*_C = 12.46 Å. However, the result from our new Eq. 24 uses D₁₅, which is right at the limit where the square root becomes imaginary, and so it probably gives a too small value, which happens to compensate for the positive deviations in Fig. 9. The conventional result involves even more accidental cancellations, as shown in Fig. 9.

We next estimate the full chain length L_C, using

(26)

where D₂ was estimated by D₃, and where the extra D₁₅ term estimates the extra contribution from the n = 16 terminal methyl end. These conventions are identical to those used two paragraphs above to obtain L_C = 13.99 Å. The value obtained using Eq. 26 and the direct simulation results for D_n is L_C = 14.01 Å. Using Eq. 24 and the simulated S_n yields L_C = 14.0 Å, and using Eq. 1 yields L_C = 14.1 Å. It therefore appears that, for these simulations, either the old or the new formula gives good estimates of L_C.

It is also worth considering such quantities for different simulations. Berger et al. (1997) reported relative positions of C₄, C₉, and C₁₄ for three different simulation runs of fluid-phase DPPC. They also reported the order parameters for each run. Our test of Eq. 24 for their best run (no. 2) is presented in Table 1. Feller et al. (1997) reported similar results, which are also used to test Eq. 24 in Table 1. Comparison values are also shown for the present simulation. In the plateau region, n = 4-9, the new Eq. 24 gives superior values for the travel when compared to the conventional Eq. 1. In the extended region, n = 4-14, the conventional Eq. 1 does better, but the relative error in Eq. 24 is less than 3%.

	HYDROCARBON THICKNESS

TOP ABSTRACT INTRODUCTION METHYLENE TRAVEL AVERAGE CHAIN LENGTH HYDROCARBON THICKNESS AREA PER MOLECULE DISCUSSION REFERENCES

To obtain a better understanding of chain packing, it is especially interesting to compare the average length of the hydrocarbon chains L_C with half the thickness of the hydrocarbon chain region D_C. The latter quantity is defined by

(27)

where A is the average area per molecule, which is 61.8 Å² in this simulation, and V_C is the volume of both hydrocarbon chains, which has been determined to be 862 Å² for this simulation, using the procedure of Petrache et al. (1997). Using Eq. 27 then gives D_C = 13.95 Å. The consistency of this procedure is indicated in Fig. 10, which shows where D_C falls on a profile of the hydrocarbon probability distribution function:

(28)

where V_n = V_CH2 for n = 2-15 and V₁₆ = V_CH3. (The flatness of the hydrocarbon probability distribution in the region 8 Å < z < 8 Å, where the headgroups and water are absent, is the criterion used to obtain the ratio V_CH3/V_CH2, and then V_CH2 comes from the simulated density of hydrocarbon in this region (Petrache et al., 1997).) Fig. 10 shows where z₁₆ and z₂ fall for the upper monolayer, as well as the probability distributions for these carbons.

View larger version (24K): [in this window] [in a new window]   FIGURE 10   Probability distribution of hydrocarbon chains, using VCH2 = 26.94 Å3 and VCH3 = 53.80 Å3. Thick solid line: Contribution from both monolayers. Dashed line: Upper monolayer only. Crosses: Average positions z2 and z16. Thin solid lines: Probability distributions for n = 2 and n = 16 (terminal methyl) in the upper monolayer. Dotted line: Headgroup (including carbonyls, glycerol, and phosphatidylcholine) probability distribution in the lower monolayer.

A major and unexpected result of this simulation is that L_C determined in the previous section is numerically very close to D_C. Although this result has often been implicitly assumed, it is not a priori correct, as emphasized by Nagle (1993). This point is illustrated in Fig. 11, which shows three caricatures, each of which packs four chains, two in each monolayer. Because of this simplicity, the average chain length in these caricatures is just L_C = (L₁ + L₂)/2. The distinguishing feature of model I is early chain termination of chain 1, so L₁ < D_C, which makes L_C less than D_C. Model I has decreasing order with increasing n because the terminal methyl end of chain 2 is more disordered than the carbonyl end, whereas chain 1 is equally disordered along the chain. The distinguishing feature of model II is interdigitation across the midplane, which makes L_C greater than D_C, and which results in increasing order with increasing n. Because the experimental order parameter decreases with increasing n, model I is clearly superior to model II. However, for this simulation, which has L_C = D_C, model I must be wrong. This has led us to propose model III in Fig. 11. Model III has both early chain terminations and interdigitation, the order parameter decreases with increasing n, and L_C = D_C, so model III is the best model for characterizing these simulations.

View larger version (20K): [in this window] [in a new window]   FIGURE 11   Three simplified caricatures, labeled I, II, and III, of packing geometry. For simplicity, only four chains are shown for each model, but more realistic distribution functions can easily be envisioned. The four numbered boxes represent the volumes occupied by each of the four chains. Ignoring upturns, the chains start from the headgroup end at z = ±DC and end in the z = 0 region. The order parameter is larger when the vertical cross section of the box is smaller. Note that these caricatures are not meant to designate sn-1 versus sn-2 chains.

	AREA PER MOLECULE

TOP ABSTRACT INTRODUCTION METHYLENE TRAVEL AVERAGE CHAIN LENGTH HYDROCARBON THICKNESS AREA PER MOLECULE DISCUSSION REFERENCES

For the present simulation, A = 61.8 Å² is obtained simply by dividing the total area of either monolayer by the number of lipid molecules in that layer. The issue is whether this value of A can be obtained using the NMR order parameters. For this simulation, which conforms to chain packing model III in Fig. 11, this is easily accomplished using Eq. 27 and the preceding estimates of L_C for D_C; this gives the satisfactory result A = 862 ų/14.0 Å = 61.6 Å².

Because model III may not apply to all bilayers, it is also of interest to evaluate a different way to obtain A that uses only the low carbon numbers in the so-called plateau region. The usual formula for this is

(29)

where V_CH2 = 26.94 Å³ is the volume per methylene group in this simulation, and the factor of 4 accounts for two chains per lipid and for a factor of 2 in D_n, which, as defined, is really twice the travel per methylene. Fig. 12 shows the results from Eq. 29 for A_n versus n. These values are all too low in the plateau region n = 3-8.

View larger version (13K): [in this window] [in a new window]   FIGURE 12   Estimate of An, using Eq. 29 () and Eq. 31 (). The horizontal line indicates the actual simulated area A = 61.8 Å.

There is a basic mathematical flaw in Eq. 29, which is corrected by taking the average of the instantaneous area, A_n = 4V_CH2/D_n; this gives (Brown, 1996)

(30)

The relation

(31)

shows that 1/D_n is larger than 1/D_n because D_n² is larger than D_n². The order parameters are used to obtain D_n (Eq. 24) and D_n² (Eq. 13). Fig. 12 shows that A_n obtained using Eq. 31 in the plateau region is in reasonable agreement with the actual A if one averages A_n over a plateau region defined to be n = 3-8.

Although the average over n = 3-8 gives the correct A, it is important to understand why there is a slope in A_n. Fig. 3 shows the average positions z_n, and Fig. 10 shows that there is a significant fraction of the headgroups at values of z_n for n = 3-5. Because these headgroups take up volume and area that is not accounted for by Eq. 30, this accounts for the smaller values of A_n in Fig. 12. The probability of headgroups decreases sharply at z_n for larger values of n, but Fig. 10 shows that the probability of terminal methyls increases. These early chain terminations mean that there is more total area that is shared by fewer chains, which is also not taken into consideration in Eq. 30; this accounts for the larger values of A_n in Fig. 12. The chain termination consideration was known previously (Nagle, 1993), and this motivated considering only the plateau region. This simulation shows that the headgroups must also be considered and that there is hardly any region that can be said to be free of both artifacts. Although one can obviously define a plateau region that works for this simulation, it is not clear if this will be universal for other simulations or for real lipid bilayers.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHYLENE TRAVEL AVERAGE CHAIN LENGTH HYDROCARBON THICKNESS AREA PER MOLECULE DISCUSSION REFERENCES

The main result in this paper has been the development of Eq. 24, which gives a new approximation for obtaining methylenic travel D_n from NMR order parameters S_n. As shown in Fig. 1, this formula fits the results of this simulation better than older approximations. Of course, the derivation of Eq. 24 used detailed results of the probability distribution shown in Fig. 8, and so Eq. 24 can only be expected to be as good as the simulation. There are results that show that the  angle for the chains did not come to equilibrium in the simulation (Fig. 5), although this seems to make only a minor difference in the tests in Figs. 6 and 7. Furthermore, the results from the different simulation of Berger et al. (1997) are in quantitatively good agreement, and results for chain fragments from three simulations (Table 1) are in good agreement. Finally, the use of Eq. 24 to predict the chain length is in excellent agreement with the actual chain length L_C. Nevertheless, Eq. 24 is a sufficiently radical departure from standard practice that it should be tested on new and independent simulations of bilayers.

A major surprise is the simulation result that half the mean thickness of the hydrocarbon chain region D_C is numerically close to the average chain length L_C. Even though this has been an implicit assumption in NMR studies (Schindler and Seelig, 1975; Thurmond et al., 1991), it is not a priori necessary. In particular, the occurrence of early chain terminations (model I in Fig. 11) would, by itself, require L_C to be smaller than D_C. However, the effect of early chain terminations appears to be canceled by the effect of interdigitation, as illustrated by model III in Fig. 11.

The result D_C = L_C leads to a simple way of estimating area per molecule A using Eq. 27. This reproduces the actual A in this simulation. This would be a very nice result, except that it leads to a disagreement. This procedure gave A_DPPC = 71.7 Å² for DPPC at 50°C (Thurmond et al., 1991). Of course, the NMR study used the older Eq. 1, but we have shown in the third section of this paper that, for this simulation, the significant deviations in this formula for individual carbons n (see Fig. 9) average out, so that the total chain length L_C is quite well approximated by either Eq. 1 or Eq. 24. Therefore, the results in this paper support the method used by Thurmond et al. (1991). However, their result for A_DPPC is 14% larger than the value A_DPPC = 62.9 Å² obtained from x-ray studies (Nagle et al., 1996).

We have examined this disagreement further for DMPC. Two recent studies have agreed that A_DMPC = 59.6 ± 0.5 Å² for DMPC at T = 30°C. Our study (Petrache et al., 1998) used the same x-ray method as used for DPPC. The other study (Koenig et al., 1997) combined x-ray and NMR results, but only changes in L_C were obtained from NMR order parameters because of the concern that no formula was adequate to give absolute numbers. Because these two studies used quite different methods and assumptions, we believe that the common value of A_DMPC that was obtained is a reliable benchmark for testing present NMR formulae.

We used the order parameters obtained by Koenig et al. (1997) at T = 30°C in Eqs. 24 to estimate D_n, which were then added up, using Eq. 26, to obtain L_C = 11.95 Å. We then assumed that D_C = L_C in Eq. 27. We also used V_C = 782 Å³, which was obtained by subtracting the headgroup volume V_H = 319 Å³ (Sun et al., 1994) from the lipid volume V_L = 1101 Å³ (Nagle and Wilkinson, 1978). The result obtained with Eq. 27 is A_DMPC = 65.4 Å², 10% higher than obtained previously. Using the conventional Eq. 1 instead of Eq. 24 reduces this marginally to A_DMPC = 64.3 Å². Using the plateau method described in the previous section and averaging A_n over n = 3-8 gives A_DMPC = 65.2 Å². Therefore, for DMPC, the methods that work so well for this simulation give values for A_DMPC that are too large.

The resolution of this disagreement may involve several issues. The first and most obvious is whether this simulation is misleading. Clearly, other simulations should be analyzed with NMR interpretation in mind. This involves testing Eq. 24, which also has implications for changes in chain length used by Koenig et al. (1997). It also involves testing L_C = D_C and Eq. 27. Of course, it is possible that all simulations will give the same but not the correct answer to the issue of equality of L_C and D_C because chain packing, as contrasted with chain conformational order, does not equilibrate in the nanosecond time range. In this regard NMR studies indicate slower, collective motions (Nevzorov et al., 1998). It may be noted that the simulation studied in this paper was a particularly long one, although a new hybrid Monte Carlo/molecular dynamics method may, in the future, overcome some equilibration barriers (Clark et al., 1999). Perhaps one might also question whether, because of subtle effects of relative time scales (Brown, 1996), NMR may not measure the same instantaneous order parameters that are obtained straightforwardly from simulations. Measured NMR order parameters that are only ~10% too low would also account for the disagreement.

We have not yet achieved the desired goal of being able to advocate an analysis method that can use NMR order parameters to quantitate bilayer structure. Nevertheless, flaws in previous