Diffraction-based density restraints for membrane and membrane-peptide molecular dynamics simulations

Abstract

We have recently shown that current molecular dynamics (MD) atomic force fields are not yet able to produce lipid bilayer structures that agree with experimentally-determined structures within experimental errors. Because of the many advantages offered by experimentally validated simulations, we have developed a novel restraint method for membrane MD simulations that uses experimental diffraction data. The restraints, introduced into the MD force field, act upon specified groups of atoms to restrain their mean positions and widths to values determined experimentally. The method was first tested using a simple liquid argon system, and then applied to a neat dioleoylphosphatidylcholine (DOPC) bilayer at 66% relative humidity and to the same bilayer containing the peptide melittin. Application of experiment-based restraints to the transbilayer double-bond and water distributions of neat DOPC bilayers led to distributions that agreed with the experimental values. Based upon the experimental structure, the restraints improved the simulated structure in some regions while introducing larger differences in others, as might be expected from imperfect force fields. For the DOPC-melittin system, the experimental transbilayer distribution of melittin was used as a restraint. The addition of the peptide caused perturbations of the simulated bilayer structure, but which were larger than observed experimentally. The melittin distribution of the simulation could be fit accurately to a Gaussian with parameters close to the observed ones, indicating that the restraints can be used to produce an ensemble of membrane-bound peptide conformations that are consistent with experiments. Such ensembles pave the way for understanding peptide-bilayer interactions at the atomic level.

FIGURE 1

Images of the DOPC bilayer plus melittin (MLT) system taken from a frame of the molecular dynamics simulation of the system. The side view is along an axis parallel to the membrane plane. The top view is along the bilayer normal. To avoid problems that might arise from the use of periodic boundary conditions, the MLT peptides in α-helical conformation were initially placed in a criss-cross pattern as shown in the top view to minimize the distance and interaction between the peptides in the two leaflets.

FIGURE 2

Images of the “box” of liquid argon before and after the application of the restraints (Table 1). Four restraint groups were randomly chosen throughout the box with the argon atoms colored according to group. Each group had the same distribution target width, but differed in the distribution target mean positions. After 40 ps of dynamics with the restraints, the four groups effectively separate into distinct distributions according to the restraint parameters.

FIGURE 3

Evolution of the distribution parameters for the four restraint groups in the liquid argon example (Fig. 2). The evolution of the group mean positions and widths are shown in panels A and B, respectively. The dashed lines indicate the target values, and the solid lines, the simulation values (colored according to the scheme indicated in Fig. 2). As the center of the box was initially at the center of the simulation cell, the groups with mean restraint values closer to the center reach their respective targets more quickly. However, all of the groups rapidly reach their target restraint values within 40 ps.

FIGURE 4

Evolution of the double-bond and water group distribution parameters for the pure DOPC bilayer system. The mean positions for the double-bond and water groups are shown in panels A and B, respectively, with the upper- and lower-leaflet values differentiated by the solid blue and red lines. The absolute value of the lower-leaflet results is shown to place both data sets on the same scale. Similarly, the distribution widths for the double-bond and water groups are shown in panels C and D, respectively. The target values are indicated by the dashed lines. During the course of the simulation, the force constants on the mean position and width terms of the restraints were increased until the target values were reached. Table 1 gives the final force constant values used. During the first 6 ns, the force constants were adjusted to bring the distribution parameters to their target values, as indicated by the regions separated by vertical bars in panel B (though these regions are only shown in B, they are the same for all of the panels). After 6 ns, the target distribution values were reached and the remaining 5 ns were used in the analysis of the trajectory.

FIGURE 5

Evolution of headgroup component group mean positions and deuterium order parameters for lipid chain carbons. (A) Mean positions for the choline (dark blue), phosphate (red), glycerol (green), and carbonyl (light blue) groups are shown. (B) Evolution of the 〈S_CD〉 for the C6 (blue) and C13 (red) carbons on the lipid tails. The average values from the last 5 ns are indicated by the dashed lines in both panels. Stable mean position and 〈S_CD〉 values are reached after 6 ns of simulation time.

FIGURE 6

Eight-order neutron and x-ray Fourier reconstructed scattering-length density profiles and structure factors for the experimental (green), unrestrained (blue), and restrained (red) pure DOPC system. (A and B) Neutron results. (C and D) X-ray results. The restrained neutron scattering-length density profile shows better agreement with experiment compared to the unrestrained system in both the density profile and the structure factors. The x-ray profile, however, shows more deviations from the experimental profile, due in part to the larger difference in the d-spacing of the restrained system compared to the experimental system (3.9 Å) compared to the unrestrained system with the experimental system (1.3 Å). Although the x-ray structure factors for both the restrained and unrestrained systems are similar, the differences in the density profiles can also be attributed in part to the difference in their d-spacings.

FIGURE 7

Evolution of the d-spacing and area/lipid for the pure DOPC simulation. With the density restraints applied to the double-bond and water groups, the d-spacing (black solid line) and area/lipid (gray solid line) initially rapidly change from the starting values, then slowly reach stable values at 6 ns as the rest of the bilayer structure equilibrates. The average values of the last 5 ns of the simulation are indicated by dashed lines.

FIGURE 8

Comparison of NVT simulation x-ray and neutron scattering-length density profiles with the experimental profiles. The x-ray profiles are shown in panel A, and the neutron profiles are shown in panel B, with the solid lines indicating the experimental results and the dashed lines indicating the simulation results. Even when the experimental membrane surface area is used in the simulation, differences between the simulation and experimental profiles are still present.

FIGURE 9

Water distribution calculated from the NVT simulation with the corresponding Gaussian fit. The direct simulation results are shown by the black solid line, and the Gaussian fit is shown by the gray dashed line. Based upon the Gaussian fit parameters, given in Table 2, the simulation water distribution is shifted and wider than the experimentally determined distribution.

FIGURE 10

Evolution of the distribution mean positions and widths for the two MLT restraint groups. The distribution means are shown in panel A, and the distribution widths in panel B. The upper- and lower-leaflet results are differentiated by the line colors and the dashed lines indicate the target values. The restraint-force constants were not varied over the course of the simulation. After 6 ns, the target distribution values were reached. As with the other simulations reported here, the distribution widths required more time to reach the target values compared to the mean positions.

FIGURE 11

Evolution of the individual MLT distribution parameters in each leaflet. In the upper leaflet ((A) distribution mean and (B) distribution width), the two MLT peptides behave similarly (shown separately in blue and red solid lines). In the lower leaflet, however ((C) distribution mean and (D) distribution width), the two MLTs appear to sample slightly different distributions, indicating that the restraints do not simply fix the proteins in the bilayer, but, rather, allow them to explore configuration space within the range determined by the experimentally-derived density distributions. The average values of the parameters are indicated by the dashed lines.

FIGURE 12

Five-order Fourier reconstructions of scattering-length density profiles. (A) Profiles for the pure DOPC bilayer without restraints (black) and the DOPC + MLT system with restraints (gray). (B) Density profiles for the DOPC + MLT system including the MLT contributions (gray) and excluding the MLT contributions (black). The MLT distribution is shown by the black dashed line. Comparison of the DOPC + MLT density with the pure bilayer density in A shows that addition of the peptide to the bilayer perturbs the pure bilayer structure, which was also observed experimentally (15). To determine the MLT distribution from the simulation, the difference structure factors were computed from the DOPC + MLT system with and without the MLT contributions.

FIGURE 13

Analysis of the MLT distribution. (A) Five-order Fourier-reconstructed distribution (solid line) along with a Gaussian fit (dashed line), which accurately describes the data. (B) Density profiles obtained directly from the simulation (solid lines) are shown along with the Gaussian fit obtained in panel A (dashed lines). There are slight differences in the two views of the MLT distribution, but both are qualitatively the same. Finally (C), the individual MLT distributions (solid lines) are shown with the total MLT distribution (dashed lines). Although the MLT peptides in the lower leaflet appear to sample slightly different distributions, overall, the individual and total distributions are very similar.

FIGURE 14

Comparison of simulation x-ray and neutron scattering-length density profiles reconstructed at the simulation and experimental d-spacing values, with experimental profiles. (A) X-ray simulation profile reconstructed at the simulation d-spacing (45.2 Å, dark gray solid line) shows more differences with the same reconstruction at the experimental d-spacing (49.1, gray solid line), compared to the corresponding neutron reconstructions from the simulation data in panel B (simulation d-spacing, dark gray solid line; experimental d-spacing, gray solid line). In both panels, the experimental profile is shown by the dashed line for reference. These profiles illustrate the increased sensitivity of x-ray scattering-length densities to d-spacing, compared to neutron scattering-length densities on the per-lipid scale.