A. Translational Diffusion
The uncorrected MSD versus time from 0 to 100ps is non-linear (, bottom) with a higher slope compared to the linear behavior at t > 1ns (, top). This is consistent with the fast diffusion region over short distances and times observed in incoherent quasielastic neutron scattering (QENS).24,28
This motion reflects rapid fluctuations of the lipids as they “rattle in a cage” of their nearest neighbors,24
and is described by the short-time diffusion constant
. It is distinct from long-time lateral diffusion, which entails hopping between
cages and is described by D
Figure 1 The lateral mean squared displacement of the lipids, l(t)−l(0)2, for 72 PME (an average of three independent simulations) and one 288 PME simulation for all 50ns (top) and for the interval 0-100ps (bottom). Standard errors are (more ...)
lists the values of
from the MD simulations. Results for both the PME simulations are in quite good agreement with experiment and with the simulated values of Lindahl and Edholm7
and Patra et al.8
The long-time diffusion constants are an order of magnitude smaller than
. These will be considered in detail for the PME systems after a brief discussion of the results from the force-shift simulation.
for the 72 FSW are substantially smaller than those from the PME simulations ( and ). The reason for this is clear from , which shows the final (50 ns) snapshot of each system. The PME bilayers (left and middle panels) are in the state of disorder expected for the liquid crystalline phase. In contrast, the aliphatic chains in the 72 FSW system are mostly fully extended and interdigitated (, right). The observed restricted translation would be expected for such a system. Nucleation of this configuration began at approximately 10ns, and condensation proceeded rapidly. It is a consequence of unbalanced electrostatic forces and the constant area constraint. Had the system been simulated at constant surface tension, the area would have most likely rapidly decreased thereby preserving a bilayer structure, albeit a highly ordered one. Previous NPT simulations with a larger electrostatic cutoff (twin range of 10/18 Å) result in a DPPC liquid crystalline bilayer29
also point to the short electrostatic cutoff as the underlying cause for the observed interdigitation. The present structure is similar to the minor domain (m
) of the DPPC ripple phase as determined from experiment by Nagle and co-workers30
and later supported via simulation by de Vries et al.31
However, it is an unphysical artifact under the present conditions (the ripple phase for DPPC is stable between of 308 to 315K32
), and further analysis is not included.
(color online) Final snapshot of one 72 PME (left), 288 PME (middle), and 72 FSW (right) simulation. The ends of the aliphatic chains are shown in large yellow spheres.
Turning to the long-time translational diffusion of the 288 lipid PME system, D
/s. This is reasonably close to results from fluorescence recovery25
and magic angle scattering NMR experiments33
(). It is also close to the diffusion constants evaluated from two GROMACS-based simulations: 1.2×10−7
/s by Lindahl and Edholm7
(with twin range spherical cutoffs, 64 lipids), and 1.3×10−7
/s by Patra et al.8
(with PME, 128 lipids). The latter authors also showed an over ten-fold dependence of D
on the spherical cutoff value for the 128 lipid system, but did not examine other system sizes.
In contrast to
, long-time lateral diffusion depends on the system size, i.e., D
/s for 72 PME, over three times larger than 288 PME. An analysis of this substantial finite size effect begins with the COM trajectories of individual lipids. shows one such trajectory from a 72 PME system. Fluctuations of the COM in relatively localized regions are punctuated by transitions to other regions, in what appears to be a combination of rattling in a cage and jumping between cages. The cluster analysis described in the Methods yields the 5 clusters for this lipid (depicted as circles in ). The average distance between clusters, d
, for all lipids is 9.05 and 8.23 Å for 72 PME and 288 PME, respectively. Given that the surface area per lipid is 64 Å2
, these values of d
and the general features of suggest that a jump model is suitable for the present analysis.
(color online) The lateral trace (in 10ps intervals) of the center of mass (COM) of a lipid in one 72 PME simulation, with the five clusters with radius of 7 Å.
plots the jump times for one of the 72 PME systems (top) and 288 PME (bottom). The 72 PME system contains periods of relative quiet and bursts of highly correlated transitions; some sets of transitions appear to be coupled between leaflets. The 288 PME system does not exhibit such obvious trends, though some degree of correlated motion cannot be ruled out. Diffusion constants estimated from the jump model,
, show a substantial system size effect ().
also differs significantly from D
for 72 PME, implying that the jump rate in this system is not well described by an independent Poisson distribution.
is slightly lower than D
for 288 PME, but the difference is not statistically significant.
Jump times (defined by the clustering algorithm and denoted by vertical tick marks) of all of the lipids in 72 PME (top) and 288 PME (bottom).
compares the probability distribution, P(Nj), for number of lipid jumps (Nj) per leaflet in two ns intervals observed in the simulations and a Poisson distribution with the same mean. The average number of jumps per 2ns block, λ, is larger for 288 PME as expected from the system size. However, the increase in λ is only two-fold because the diffusion is slower than 72 PME. The probability distribution for 72 PME (top panel of ) appears to be bimodal. This is consistent with the correlated bursts of activity shown in , and qualitatively different from a Poisson distribution. P(Nj) for the 288 PME is closer to a Poisson distribution than 72 PME (, bottom panel) although is somewhat broader, implying some correlation in lipid jumps even for this system.
The average probability of the number of jumps in a 2ns block (λ) ignoring the first and last block for 72 PME (top) and 288 PME (bottom), and the Poisson distribution for these values of λ (solid lines in each panel).
Correlation, or concertedness, in the jumps was analyzed by the approach of Brown et al.27
described in the Section II-B. lists the relevant statistics for the 100ps window size. The individual 72 PME systems show substantial evidence of concertedness in the first and second shells (p
≤ 0.1 in all three cases), and run1 even shows potential correlation to the third shell. Similar results were observed for the 50 and 200ps time windows, but no concertedness is evident for the 500ps time window. Pooling the 100ps time window results of the three simulations yields p
for the first and second shells, and p
= 0.024 for the third shell. Hence, there is little doubt that concertedness extends to the second shell and it likely extends to the third shell for systems of 72 lipids.
Table 2 Statistics for concerted motion: N, the number of diffusional jumps for each trajectory based on a cluster analysis; Nf, the number of jumps following another jump within 100ps; Nif, the number of jumps that would be expected from independent Poisson (more ...)
Results for 288 PME are somewhat less clear, partly because of sample size. There is no strong evidence for concertedness with the first shell (p = 0.31), but correlation with the second shell is high (p = 0.005). Correlation with the third and fourth shells is negligible (data not shown). Similar results are obtained with the 50 and 200ps time windows; the 500ps window yields no correlation for any of the shells. It is not clear at this juncture whether the lack of correlation for the first shell arises from incomplete sampling, limitations in the jump model, or, more deeply, to the mechanism of lipid diffusion. Nevertheless, the presence of correlation with the second shell is the central result of the study. If diffusion exhibits a natural correlation with the second shell (a length of 16 Å), it becomes likely that the effects of a transition of a lipid in a system of 36 lipids/leaflet will propagate across the 50 Å/side box length, and thereby artificially increase the diffusion constant. Because the 288 PME system is nearly 200 Å/side, such artifacts are unlikely.
The finite size effect reported here is opposite of that in pure bulk simulations, where the self-diffusion constant increases with system size.4,5
Yeh and Hummer5
determined that the effect in homogeneous systems is hydrodynamic in origin. In bilayers, the present finite size effect is shorter range and predominantly collisional. It is possible that a longer-range hydrodynamic effect will be observed in simulations of larger bilayers.
Recent simulations of a 1-stearoyl-2-oleoyl-phosphatidylethanolamine (SOPE) lipid bilayer of 72 lipids showed hydrogen bond networks that span the simulation box size.34
In contrast to the present DPPC simulations, the calculated D
is lower than experiment by a factor of about three. While the time spent in the clusters is small, such the hydrogen bond networks might be expected to retard diffusion. Simulations SOPE bilayer with larger sizes are needed to determine whether such large hydrogen bond networks are partly finite artifacts and are the cause of the large underestimate of D
Lastly, if the COM displacement of each monolayer is removed from the lipid COM displacement, the lateral diffusion constant (
) is reduced by more than factor of 4 for 72 PME, but only 6% for 288 PME (). The origin of this finite size is two-fold. First, as the size of the bilayer increases, the effect of lipid diffusion on the COM of each monolayer decreases; i.e.,
. Second, correlated motion leads to substantial unidirectional displacements for 72 PME. Removal of the monolayer COM removes some of these, further reduces
, and the errors partially cancel. In general, monolayer COM displacement should not be removed when evaluating diffusion constants of lipids in bilayers.
B. Finite size effects on bilayer structure and lipid reorientation
The electron density of the lipid bilayer and deuterium order parameters of the lipids are equal for our 72 and 288 PME simulations. These results agree with de Vries et al.3
, where structural properties converged at a simulation size of 72 lipids (36 per leaflet). Although these structural properties remain unchanged, the enhanced lipid diffusion calculated for the 72 PME may influence radial and angular distributions of the lipid charged atoms. Therefore, the two-dimensional radial distribution function, g2D
), for the P and N atoms in DPPC for 72 and 288 PME is now considered. The g2D
) averages out the bilayer normal axis and is a distribution function per leaflet. plots g2D
) for P-N and N-N and demonstrates that the distributions are identical for both system sizes. The distribution of the relative orientation of the P-N dipole is also a measure of changes in angular distributions. shows the angular probability of all lipids within 7 Å of the reference lipid (averaged over all lipids). A cutoff value of 7 Å was chosen because it resulted in a distribution with the largest peak at 180°, i.e., neighboring dipoles are in opposite direction. The difference between 72 PME and 288 PME is negligible. These results demonstrate that even though the lateral diffusion shows system size dependence the head group radial and angular distribution functions do not.
The 2-dimensional radial distribution functions, g2D(r), for one leaflet for P-P (top) and N-N (bottom).
The probability of finding a neighbor with an angle θ defined by the angle formed between its P-N vector and those neighboring lipids less than 7 Å of separation between the phosphorous atoms.
The hopping rate in the 72 PME systems is approximately 0.1ns−1
(). Consequently, rotational motions of the lipid, which take place on the 1 to 10 ns time scale,35,36
could be perturbed by the finite size effect observed for lateral diffusion. To investigate this possibility, three vectors were considered: the principal axis (PA), phosphorous-nitrogen (P-N) vector, and phosphorous-hydrogen (P-H) vector with the neighboring hydrogen on glycerol. The P-N vector is a measure of head group motion, PA is associated with wobble,35
and P-H with axial diffusion.36
plots the reorientational correlation functions for these three vectors for 72 and 288 systems. There is a slight deviation for the C2
) of the PA, but the three runs of 72 PME bracket the 288 PME result (only the average is shown). The fastest relaxing P-N vector shows negligible differences between the small and large systems.
Figure 8 (color online) Second rank reorientational correlation functions of the principal axis (PA), P-N, and the vector between the phosphate and adjacent glycerol hydrogen (P-H) (top), and C-H vectors of aliphatic carbons C2, C9, and C15 (bottom) for 72 PME (more ...)
Next, the fast torsional motions of the carbon-hydrogen vectors on the aliphatic chain are compared for the two systems. The reorientational correlation functions for the C2, C9, and C15 positions (C2 is closest to the headgroup) are shown in the bottom panel of . Similar to the above mentioned correlation functions, the decays are nearly identical for the 72 PME and 288 PME, with the expected decrease of relaxation times with increasing distance from the head group. The invariance of C2(t) to system sizes greater than or equal to 72 lipids allows for force field parameterization on torsional changes to the potential function based on results of a smaller system of 72 lipids.