Molecular Simulations of Mixed Lipid Bilayers with Sphingomyelin

Apr 27, 2017 - ... and ‡Biophysics Program, University of Maryland, College Park, .... Solid-State Deuterium NMR Spectroscopy of Membranes. Trivikram ...
0 downloads 0 Views 6MB Size
Article pubs.acs.org/JPCB

Molecular Simulations of Mixed Lipid Bilayers with Sphingomyelin, Glycerophospholipids, and Cholesterol Indrani Bera† and Jeffery B. Klauda*,†,‡ †

Department of Chemical and Biomolecular Engineering and ‡Biophysics Program, University of Maryland, College Park, Maryland 20742, United States S Supporting Information *

ABSTRACT: The highly diversified composition of lipid bilayers across living cells is crucial for many biological processes. Lipid bilayers mainly consist of phosphatidylcholines (PC), phosphatidylethanolamines (PE), sphingomyelin (SM), and cholesterol, with eukaryotic membranes containing high percentage of sphingomyelin and cholesterol. In this study, we have modeled bilayers with different concentration of PC, PE, and SM to understand the changes in bilayer properties with varied SM concentrations. In addition, membrane models with 33% cholesterol have been simulated to understand the influence of cholesterol. To quantitatively access the structure and dynamics of membranes, deuterium order parameters (SCD), mass density profiles, lipid relaxation times, clustering analysis, and radial distribution functions are calculated. The SCDs compare favorably with past NMR experiments and increase with an increase in SM content. The surface area calculations showed that on addition of 50% palmitoyl-SM (PSM) surface area decreases (60.0 ± 0.6 Å2) from that of pure POPC (64.7 Å2), which is further lowered in the presence of cholesterol (44.4 ± 0.2 Å2). The lipid axial relaxation time decreases with increase in concentration of glycerophospholipids. The accuracy of these lipid membranes allows for future studies with more complex lipid mixtures containing SM to represent the diversity of lipids in natural membranes.



INTRODUCTION Lipid molecules that can assemble into bilayers about 5 nm thick at certain water-to-lipid ratios provide the basic structure for all cellular membranes. The outer cell membrane serves as a wall to protect the cell, and transmembrane proteins act as gatekeepers for passage of different solutes inside the cell. Activities of transmembrane proteins and other membraneassociated proteins are governed by the composition of the lipid bilayer. Phosphatidylcholines (PCs), phosphatidylethanolamines (PEs), sphingomyelins (SMs), and cholesterols are common lipids present in eukaryotic cell membranes.1 The abundance of these lipids varies with different cell membranes that exist in different organisms and organelles.2 SMs and cholesterols are present in the outer leaflet of the plasma membrane.3 From studies over the past few decades, it has been found that the structural organization of lipid bilayers, their phase behavior, and mobility are governed by the presence of SMs and cholesterols.4,5 At the correct thermodynamic conditions, cholesterol−lipid interactions result in phase separation into a liquid-ordered and liquid-disordered phase in model membranes of a few lipids. These phases are macroscopic in nature and can extend up to several micrometers.6 The presence of SMs and cholesterols in complex mixtures and natural membranes results in inhomogeneity in the bilayer as they tend to form local clusters and quite often results in the formation of lipid rafts.7−9 The lipid © 2017 American Chemical Society

rafts are laterally segregated nanodomains in the biological membranes containing high amounts of SMs and cholesterol,10 and they are known to influence many biological processes11,12 and are very difficult to detect due to their size and dynamic nature. These rafts play important role in sorting of membrane molecules and signal transduction.10 An increasing amount of cholesterol results in tighter packing of the lipid bilayer and also influences the function of membrane proteins.7,8 The lipid bilayers forming rafts are clinically very important as raftassociated proteins are targets for many diseases such as Alzheimer’s disease, Parkinson’s disease,13 and muscular dystrophy.14 The functions of the biological membranes largely depend on the structure and composition of bilayers. It is very difficult to resolve atomic structures of lipid membranes and membrane proteins using experimental approaches. Molecular dynamics (MD), on the other hand, is widely used as a complementary method to resolve atomic coordinates of lipids in membranes and study the membrane structure and dynamics.15,16 Initial MD studies on lipid bilayers involved only one lipid,17,18 whereas later on, many studies on mixture of lipids have been reported from simple mixtures19−22 to those that represent Received: January 12, 2017 Revised: April 26, 2017 Published: April 27, 2017 5197

DOI: 10.1021/acs.jpcb.7b00359 J. Phys. Chem. B 2017, 121, 5197−5208

Article

The Journal of Physical Chemistry B

Figure 1. Chemical structures of the lipids used in the study.

various membranes in organisms and their organelles.23−30 Several simulation studies have been done to investigate the properties lipid mixtures composed of sphingolipids and cholesterols31,32 and also cholesterol and PCs.33−37 Comparative simulations studies for SMs and PCs showed that SMs have lower surface area per lipid than the PCs.38 Thus, it is expected that a bilayer consisting of both SMs and PCs will have lower surface areas per lipid than pure PCs. The addition of cholesterol to either SMs or PCs resulted in lower surface areas than the individual ones in simulations.31 These initial simulation studies were carried out using force fields that have not been tuned and verified with updated experimental data on SM lipid structure.39,40 To develop proper membrane models, consistent with experimentally determined structural properties, accurate force field parameters are necessary. Recently, improved and more accurate all-atom force fields have been developed for PC, sphingolipids, and cholesterol.41−43 A recent study has provided a critical comparison of using different force fields for studying structure and dynamics of simple lipid bilayer (DMPC, POPC, and POPE).44 The study shows that the GROMOS54a7 united-atom force field underestimated the lipid self-diffusion and could not describe the melting characteristics of lipid bilayers, whereas CHARMM36 behaved better. Slipids (an all-atom force field) is well suited for studying membrane bilayers above the melting temperature. The study also showed Lipid14 and CHARMM36 could define well-ordered membrane structures for the gel phase and a highly cooperative transition during melting. Moreover, we have shown that the GROMOS force field for sphingolipids are likely too dense in comparison with recent NMR experiments.42 This does not mean that united-atom models are fundamentally flawed, but rather demonstrates a need to reoptimize parameters with updated experiments. Only a recent study using CHARMM36 force field has been

performed on ternary mixtures of PCs, sphingomyelin, and cholesterol at concentrations depicting liquid-ordered, liquiddisordered, and coexistence of the two phases.45 In this study, we have generated binary mixtures of lipid bilayer models consisting of PCs and PEs with SMs and ternary mixtures containing cholesterol along with these lipids using the most recent CHARMM36 lipid force fields.



METHODS Bilayer Setup and Simulation Protocol. Since concentrations of different lipids result in different structural properties of the bilayer and consequently effect the biological function, modeling of bilayers with SM and cholesterol are necessary to understand the effect of the presence of these lipids on the PC- and PE-rich membranes. In this work, all the PE and PC lipids were monounsaturated with unequal sn-1 and sn-2 chains. We have used the 16:0 and 18:1 fatty acids as hydrophobic tails in our studies along with above-mentioned head groups. Also, we have built bilayers with different concentrations of PC and SM to understand the effects of changes in concentrations to the structure of the bilayer. We generated binary system models with 1-palmitoyl-2-oleoyl-snglycero-3-PC (POPC), 1-palmitoyl-2-oleoyl-sn-glycero-3-PE (POPE), and N-palmitoylsphingomyelin (PSM) in different concentrations and ternary systems using cholesterol and phospholipids. The structures of the lipids are shown in Figure 1, and the details of different systems are provided in Table S1. The bilayer models were generated using the CHARMMGUI Membrane Builder.46−48 All the bilayers systems were built in rectangular box containing a total of 100 lipids (50 per leaflet) for binary systems and 300 lipids (150 per leaflet) for ternary systems with a hydration number of 30 water molecules per lipid. The CHARMM TIP3P water model49,50 was used 5198

DOI: 10.1021/acs.jpcb.7b00359 J. Phys. Chem. B 2017, 121, 5197−5208

Article

The Journal of Physical Chemistry B along with the CHARMM36 lipid force field.41−43 The CHARMM TIP3P water model is a modified version of original TIP3P water model49 in which hydrogens are provided with an additional van der Waals term to prevent close contacts with other charged atoms and avoid unrealistically large electrostatic energies.50 MD simulations were carried out using GROMACS 5.0.451,52 using developed CHARMM-GUI equilibration scripts for GROMACS. The systems were energy minimized using 1000 steps of steepest descent algorithm followed by six step equilibration runs, of which first two steps composed of constant number of atoms, volume, and temperature (NVT) ensemble and the remaining steps in the NP(pressure)T ensemble.53 At least 300 ns was used for production runs in the NPT ensemble at 321.15 K and 1 bar with three replicas per system. The Nosé−Hoover thermostat was used for temperature coupling.53 For pressure coupling, the Parrinello−Rahman method54 was used with semi-isotropic scaling at pressure of 1 bar. The simulations were carried out using periodic boundary conditions (PBC), and particle mesh Ewald (PME) was used for electrostatic interactions.55 A forceswitch cutoff approach was used for the Lennard-Jones potential over 1.0−1.2 nm. The LINCS algorithm was used to restrain the bond lengths.56 A time step of 2 fs was used and coordinates were saved every 2 ps. All runs were replicated thrice to avoid biased results. For each of the three runs, a different bilayer model was generated using CHARMM-GUI. Calculation of Membrane Properties. After the production runs, the structural and dynamical properties of membranes were calculated from the equilibrated bilayers. For determining whether the lipid bilayers are equilibrated, calculation of surface area per lipid (SA) is an important metric. The stability of the systems can be assessed from the SA (details provided in the Results section). The SA was calculated by the area of the system in x−y plane divided by the number of lipids in each leaflet. The component areas were calculated using a triangulation method in Qhull.57 The method uses coordinates from lipids to generate a polygon from which the surface area of the lipid is calculated. Three atoms from phospholipids and one atom from sterol (C2, C21, and C31 for glycerol lipids; C2S, C1F, and C4S for sphingo lipids and O3 for sterols) are selected as Voronoi vertices. Another important metric for prediction of structural of lipid bilayers is deuterium order parameters (SCD).58,59 In NMR experiments, SCDs are largely used to determine the order of acyls chains of lipids,60 which give insights about the phase and configuration of the lipids in the bilayer. The order parameters are expected to decrease along the length of the acyl chain. SCDs are calculated using the equation SCD =

3 1 cos2 β − 2 2

as the distance between the two leaflets of bilayer’s hydrocarbon acyl.61 Inter- and intramolecular hydrogen bonds were calculated using gmx_hbond. Lipid motions can be calculated from reorientational correlation functions and can be related to relaxation times from NMR. The relaxation times are measured by second-order reorientational correlation function, C2(t), using the equation C2(t ) = ⟨P2[μ(0) ·μ(t )]⟩

(2)

where P2 is the second Legendre polynomial and μ is a particular atom vector. The correlation function studied here is the cross-chain vector for C2−C2 on each chain and describes how much the spin of the vector at time t correlates to the spin of the vector at t0. Exponential fits were obtained to determine the time constants using the following empirical function:62 3

C2(t ) = a0 +

∑ aie−t /τ

i

i=1

(3)

a0 is estimated from the average value of plateau reached by the correlation function, a1, a2, and a3 are coefficients of three exponential decay functions, τ1 and τ2 represent the internal motions, and τ3 represents the lipid wobble motions62 based on observed time scales and more physically based motional models. The radial distribution function, g(r), gives the probability of finding a particle in the distance r from another particle. Radial distribution functions were calculated in 2D using gmx_rdf. Lateral clustering of lipids was done to visualize the effects of hydrogen bonds and other intralipid interactions. The clustering of lipids were done using python scikit-learn algotithm63 which uses density based spatial clustering64,65 successfully implemented in61 with a Dcut value of 4 Å. The visualizations were done using VMD.66 All the analyses were performed on the last 150 ns of the simulation runs.



RESULTS Surface Area per Lipid. All the simulations maintained stable bilayers and sample snapshots of the systems after the 300 ns runs are in Figure 2. The stability of the bilayers after the 300 ns simulation runs were checked by calculating average SA per lipid. Nearly all the systems reached equilibration in the very initial stage (before 50 ns) (Figures S1−S4). The POPE− PSM (1:2) run2 initially appeared to be going to gel-like phase, but the extended simulations to 500 ns (Figure S3) showed the systems to be in equilibrium (not in a gel-like phase). Therefore, all further analysis on bilayer properties were taken from the last 150 ns of the simulations. The average of SA/lipid for all lipids was obtained and is listed in Table 1. It is known that for a single-component membrane that PSM has a lower SA/lipid compared to PCs.38 Also, the SA/lipid is lower for PEs as compared to PCs which has been determined using solid-state 2H NMR.59,67 Our MD simulations show that addition of PSM to the bilayer containing POPC or POPE results in more ordered phase, which is indicated by the lowering of SA/lipid as the PSM concentration is increased (Table 1). Comparing with the 1:1 PC (or PE):PSM, the effect of increasing POPC concentration is greater for the SA/lipid compared to POPE when going from the equal molar to higher PC (or PE) concentrations. When increasing the PSM concentration, the condensing effect is similar between membranes with POPC and POPE. Our simulation study resulted in SA/lipid of 44.44 ± 0.2 Å for

(1)

where β is the angle between each of the CH bond in the lipid tail and the bilayer normal. Mass density profiles (MDPs) were calculated for individual groups using gmx_density tool along the bilayer normal for individual groups such as phosphate, choline, glycerol, carbonyl, CH, CH2, and CH3 along with total mass density. Headgroupto-headgroup distance (DHH), the overall bilayer thickness (DB), and hydrophobic thickness (2DC) were calculated from the MDPs. DHH is defined as the distance between the peaks of total mass density, DB is defined as twice the distance between the midpoint of maximum water densities, and 2DC is defined 5199

DOI: 10.1021/acs.jpcb.7b00359 J. Phys. Chem. B 2017, 121, 5197−5208

Article

The Journal of Physical Chemistry B

POPE (∼6 vs ∼4 Å2). The individual SA of the lipids (POPC and PSM) decreased more than 6 Å2 on addition of cholesterol. Deuterium Order Parameters. The order parameter is important to quantify the structural orientation and flexibility of lipids in bilayers and can be determined experimentally with NMR. The carbons near the headgroup have higher order, which decreases down the length of carbon chain since the motions of head groups are restricted and the tail regions are relatively free to move in the bilayer. The SCD parameters from our simulation studies show this trend (Figures 3 and 4). Since

Figure 2. Snapshots of different systems (A) POPC−PSM1:1, (B) POPE−PSM1:1, (C) POPC−PSM−Chol1:1:1, and (D) POPE−PSM− Chol1:1:1. POPC: red; POPE: magenta; PSM: green; cholesterol: yellow.

Table 1. Component Surface Area per Lipid and Average Surface Area of the System in Å2 in the Different Bilayer Systems ± Standard Error system POPC−PSM2:1 POPC−PSM1:1 POPC−PSM1:2 POPE−PSM2:1 POPE−PSM1:1 POPE−PSM1:2 POPC−PSM− cholesterol1:1:1 POPE−PSM− cholesterol1:1:1

POPC or POPE 63.9 62.8 62.2 57.7 56.6 57.0 55.5

± ± ± ± ± ± ±

0.3 0.3 0.5 0.2 0.6 0.5 0.3

54.2 ± 0.3

PSM 58.0 57.1 58.8 53.2 53.0 53.6 54.0

± ± ± ± ± ± ±

Chol

av SA/lipid [Å2] ± ± ± ± ± ± ±

0.5 0.7 0.1 0.5 0.4 1.7 0.1

27.9 ± 0.1

62.0 60.0 58.2 56.2 54.8 53.6 44.4

52.9 ± 0.2

27.2 ± 0.5

43.5 ± 0.3

Figure 3. SCD parameters for sn-1 chain of POPC or POPE in (A) POPC−PSM1:1 and (B) POPE−PSM1:1 compared with the respective simulations of pure POPC72 and POPE41 bilayers. Orders parameters from mixed bilayer are shown in red, and those from pure systems are shown in blue.

0.6 0.6 0.7 0.7 0.1 0.9 0.2

our previous MD simulations with pure POPC and POPE bilayers demonstrate excellent agreement with NMR experiments,41,61 we expected similar agreement with mixed lipids without cholesterol. The average order parameter for fatty acid chain of PSM in POPC:PSM1:1 is 0.20, which is in agreement with experimental calculated value of 0.197.69 Figure 4 and Figure S6 compare simulation and NMR SCD parameters for POPC-PSM bilayers.39 The agreement between experiment and MD simulations with the C36 lipid force field is excellent with a slight tendency of the simulations to be higher than experiment for most comparisons. This may suggest that the SA/lipid for systems with POPC and PSM should be slightly lower than experiment due to the higher SCD. The SCDs of both sn-1 and sn-2 (Figure 3 and Figure S5) chains are higher for the mixed bilayers systems (POPE−PSM) than those of pure POPE bilayers, which correlates (inversely proportional) with the lower the SA/lipid of mixed POPE− PSM1:1 (56.6 Å2) bilayer than that of pure POPE bilayer (59.2 Å2).70 The order parameters for fatty acid chain of PSM are higher than those of acyl chain of POPC. As expected, the SCDs are higher for the bilayer with higher concentration of PSM, correlating with lower surface area per lipid molecule (Figure 4). The SCD for sn-1 acyl chain of POPC shows the same trend; the order parameters increased with increase in concentration of PSM (Figure S6). Comparing plots the SCDs of both POPE and PSM chains in POPE−PSM systems are higher than those

POPC−PSM−Chol1:1:1, which is expected because cholesterol and SMs condense bilayers and result in a more ordered phase.68 The average decrease in SA/lipid in the case of ternary mixture with POPC is more than 15 Å2 in comparison with that of binary mixtures with PSM. Interestingly, the average decrease is lower (12 Å2) in the case of mixture with POPE lipid. The SA/lipid for the individual lipid components in the bilayer systems was also calculated, and the results are provided in Table 1. We see that the component surface areas also follow with the trend of the overall SA/lipid. The component surface area of POPC decreased with an increase in the PSM concentration. The component SA of PSM decreases with an increase in the PSM concentration from 33% to 50%. There is a decrease in component SA of PSM in the presence of POPE with respect to that with POPC. This is the result of hydrogen bonding between POPE-PSM (see below). However, the difference between the component SA of the glycerophospholipid and PSM is greater with POPC compared to 5200

DOI: 10.1021/acs.jpcb.7b00359 J. Phys. Chem. B 2017, 121, 5197−5208

Article

The Journal of Physical Chemistry B

Figure 5. SCD parameters for sn-1 chain of POPC in (A) POPC− PSM−Chol1:1:1 and (B) POPE−PSM−Chol1:1:1 compared with simulations of pure POPC72 and POPE41 systems. Orders parameters from mixed bilayer are shown in red, and those from pure systems are shown in blue.

(resulting in slightly higher order parameter) than in system with POPC. Mass Density Profiles (MDPs). The DHH, DB and DC from the MDPs for different systems calculated from the averages of the three runs are tabulated in Table 2. The results show that

Figure 4. SCD parameters for fatty acid chain of the PSM lipid (A) POPC−PSM2:1, (B) POPC−PSM1:1 and (C) POPC−PSM1:2 obtained from simulations (red) and compared with experimental values (blue).39

Table 2. Average ± Standard Error Head-to-Head Distance (DHH), Bilayer Thickness (DB), and Hydrophobic Thickness (DC) for the Systems under Study

of POPC and PSM in POPC−PSM systems likely due to PE’s tendency to interact more with PSM via hydrogen bonds (see below). The order parameters increased in the presence of cholesterol due to its condensing effect on lipid membranes (Figure 5). The order parameters of sn-1 chain of POPC in ternary mixture of POPC−PSM−Chol1:1:1 were, as expected, higher than that of pure POPC (Figure 5). The same trend was observed for ternary systems with POPE (Figure 5). The SCD parameters of sn-1 acyl chain for POPE−PSM−Chol1:1:1 system compared to that of experiment is shown in Figure S7, where the experimental SCD is higher for carbons 2 and 3 adjacent to the headgroup. However, agreement for other carbons down the chain is excellent. Recently,71 other mixed lipid membranes with cholesterol and PSM were studied using NMR. Two of the membrane bilayers studied are ternary systems consisting POPC−PSM−cholesterol2:1:1 and POPE−PSM−cholesterol2:1:1 at 50 °C. Acyl chain order at carbon atom 9′ was determined using deuterium on just this carbon and not the fully deuterated systems used in other studies, as at this position cholesterol has the largest effect. The order parameter calculated for the carbon atom 9′ of POPC−PSM−cholesterol2:1:1 (25% cholesterol) was 0.336, whereas for POPE−PSM−cholesterol2:1:1 it was 0.34. In our study, which was carried out at near similar temperature of 48 °C, for POPC−PSM−cholesterol1:1:1 the same carbon has an average order parameter value of 0.378, and for the system POPE−PSM−cholesterol1:1:1 it was 0.39. Our system was with 33% cholesterol, and thus the increase in cholesterol concentration has increased the order parameters. In the presence of POPE, effect of cholesterol was slightly less

composition POPC−PSM2:1 POPC−PSM1:1 POPC−PSM1:2 POPE−PSM2:1 POPE−PSM1:1 POPE−PSM1:2 POPC−PSM−cholesterol1:1:1 POPE−PSM−cholesterol1:1:1

DHH (Å) 39.3 39.4 40.3 40.9 43.1 41.2 45.4 46.0

± ± ± ± ± ± ± ±

1.0 0.3 0.7 1.2 0.1 0.4 0.2 0.2

DB (Å) 40.3 40.3 41.6 42.5 42.9 43.2 45.7 46.1

± ± ± ± ± ± ± ±

0.5 0.3 0.3 0.2 0.8 0.3 0.1 0.1

DC (Å) 27.3 26.2 26.6 28.3 27.5 26.7 28.3 29.5

± ± ± ± ± ± ± ±

0.3 0.3 0.4 0.2 0.2 0.3 0.1 0.3

head-to-head distance and bilayer thickness for POPE−PSM1:1 are greater than those of POPC−PSM1:1. The ternary mixtures with cholesterol have a much larger bilayer thickness (>3−6 Å) and larger head-to-head distance than those of binary systems. The overall hydrophobic thickness is not affected by changing the concentration of lipids and is not influenced much but POPE systems have slightly greater hydrophobic thickness than POPC systems. The bilayers with POPE (both binary and ternary) have bilayer thickness greater than those with POPC. The bilayer thickness for pure POPC system is 37.3 Å, and the hydrophobic thickness is 28 Å.61 Thus, the addition of PSM and cholesterol increases the bilayer thickness and also increase with increase in concentration of PSM. Previous simulations with POPC:PSM:cholesterol in 62:1:1 ratio, reported head-tohead distance as 35 Å,73 whereas DHH of pure POPC was found to be 36.4 Å. Clearly, addition of PSM and cholesterol in the bilayer, the DHH increases. The total MDPs for different 5201

DOI: 10.1021/acs.jpcb.7b00359 J. Phys. Chem. B 2017, 121, 5197−5208

Article

The Journal of Physical Chemistry B

of hydrogen bonds are reported in Table 3, and the number of hydrogen bonds vs time is shown in Figure 9.

systems are shown in Figures 6 and 7, and those of the components are shown in Figures S8 and S9. The MDP of

Figure 6. Total MDP (A) POPC−PSM systems and (B) POPE−PSM systems.

Figure 7. Total MDP for ternary systems. Figure 8. Snapshots showing intramolecular (green) and intermolecular (black) hydrogen bonds between two PSM molecules (A), two different conformations of phosphate oxygen (B), and hydrogen bond between phosphate oxygen of PSM and oxygen of cholesterol (C).

phosphate and choline (PCs) of respective lipids increase with their relative concentration. The MDP of the headgroup (phosphate and choline) of PSM is farther from the center of the bilayer when compared to that of POPC, and similar comparison exists for densities of PCs in POPE when compared to that of PCs in PSM. Hydrogen Bonds. Inter- and intramolecular hydrogens bonds calculated from our MD simulations provide detail in bilayer structural properties due to molecular interactions. Sphingomyelin and cholesterol affect the bilayer properties by forming interlipid hydrogen bonds. Calculating the number of hydrogen bonds quantifies the influence of sphingomyelin and cholesterol in the bilayer. POPC forms only intermolecular hydrogen bonds with PSM, whereas POPE and PSM are involved in both intermolecular and intramolecular hydrogen bonds. The OH group of PSM at C3 is involved in intramolecular H-bond to the phosphate O atom and the amide NH forms intermolecular H-bonds to the carbonyl O atom. The amide NH also forms an intermolecular hydrogen bond with phosphate oxygen. A snapshot of two PSM molecules showing all possible inter- and intramolecular hydrogen bonds are shown in Figure 8. Average frequencies

As expected, intramolecular PSM−PSM hydrogen bonds increase with increase in concentration of PSM. The number of intermolecular hydrogen bonds between POPE and PSM is large compared to POPC and PSM due to the presence of amine group in POPE which form the additional bonds. The average number of intramolecular hydrogens for PSM is also more in the POPE/PSM system than in the POPC/PSM system. The details of intra- and intermolecular hydrogen bonds between PSM molecules themselves are shown in Table 4 and compared to that of only PSM simulation.42 In both POPC−PSM1:1 and POPE−PSM1:1, the probabilities of intramolecular hydrogen bonds between hydroxyl group and phosphate oxygen (O−H:::O−P) is reduced compared to the pure PSM bilayer. Whereas in case intermolecular PSM bonds, the probabilities of a hydrogen bond between the amide NH of one PSM and phosphate oxygen of another PSM increases and that with carbonyl oxygen decreases for the mixed systems than the pure PSM. The number of intermolecular hydrogen bonds 5202

DOI: 10.1021/acs.jpcb.7b00359 J. Phys. Chem. B 2017, 121, 5197−5208

Article

The Journal of Physical Chemistry B

Table 3. Average Frequencies of Hydrogen Bonds for Different Systems ± Standard Error Calculated for Block Size of 10 ns composition

POPC (or PE)−PSM

POPC−PSM2:1 POPC−PSM1:1 POPC−PSM1:2 POPE−PSM2:1 POPE−PSM1:1 POPE−PSM1:2 POPC−PSM−cholesterol1:1:1 POPE−PSM−cholesterol1:1:1

0.099 0.112 0.104 0.348 0.427 0.379 0.060 0.211

± ± ± ± ± ± ± ±

0.002 0.002 0.002 0.004 0.005 0.006 0.001 0.018

PSM−PSM 0.361 0.568 0.807 0.376 0.491 0.844 0.371 0.379

± ± ± ± ± ± ± ±

0.002 0.002 0.001 0.002 0.004 0.010 0.001 0.005

POPE−POPE

0.480 ± 0.007 0.310 ± 0.003 0.180 ± 0.010 0.170 ± 0.002

total no. of H-bonds/lipid 0.460 0.680 0.912 1.204 1.228 1.403 0.431 0.756

± ± ± ± ± ± ± ±

0.003 0.003 0.002 0.008 0.007 0.015 0.001 0.187

between POPE molecules should decrease which is evident for our results, though our system was studied at 48 °C. For the ternary systems, intermolecular hydrogen bonds between individual lipids and cholesterol were calculated (Table 5). We see that PSM forms more hydrogen bonds Table 5. Average Frequencies of Hydrogen Bonds between Lipids and Cholesterol for Ternary Systems ± Standard Error Calculated for Block Size of 10 ns

Table 4. Intra- and Intermolecular H-Bond Frequencies between PSM Molecules PSM42

O−H:::O−P O−H:::OP

0.99 0.01

N−H:::OP N−H:::OC

0.02 0.28

POPC−PSM1:1 intramolecular 0.87 0.011 lipid−lipid 0.08 0.16

POPC (or PE)− cholesterol

PSM− cholesterol

POPC−PSM−cholesterol1:1:1 POPE−PSM−cholesterol1:1:1

0.040 ± 0.000 0.060 ± 0.001

0.100 ± 0.000 0.110 ± 0.001

with cholesterol than the PC/PE. On average, increasing the concentration of PSM results in more hydrogen bonds per lipid (Table 3); i.e., 33% of PSM results in half of the lipids hydrogen bonded and 67% approaches one hydrogen bond per lipid with membranes mixed with POPC. PSM membranes with POPE have the largest hydrogen bonds per lipid and the effect of increasing PSM is smaller for this quantity compared to membranes with POPC. Radial Distribution Functions (RDFs) and Lipid Clustering. The RDF measures the probability of finding an atom with respect to distance from the reference atom and is a useful tool to estimate the lateral ordering of lipid bilayers. The 2D RDFs of the phosphates were calculated for all possible lipid interactions (PC/PE−PSM, PSM−PSM and PC/PE− PC/PE) and were well converged (Table S2). The first peak appears at a most probable distance where the oxygen of phosphate groups is pointing toward each other, and the second peak with lower g(r) corresponds to the conformation when the oxygen atoms are pointing away from each other. Hydrogen bonds between phosphate oxygen and amide nitrogen are formed in these conformations. For physical insight, a snapshot showing the different conformations of PSM molecules are shown in Figure 8. The RDFs for the phosphate−phosphate of PSM for different systems with POPC−PSM are shown in Figure 9. All RDF plots show two sharp peaks for PSM−PSM interactions, depicting two different conformations of the oxygen atoms (Figure 8B). The g(r) value of the first peak is slightly less for system with 2-fold concentration of POPC than for the system same amount of PSM. RDF plots for the phosphate groups of self POPC/ POPE, POPC/POPE−PSM, and self PSM for all the systems are shown in Figure S10. The g(r) value of the first peak of POPC−POPC interaction is less compared to PSM−PSM interaction, and the first peak is formed at larger distance than PSM−PSM interaction. For systems with POPE, the second peak is not well-defined for the RDFs between POPE−POPE and POPE−PSM (Figure

Figure 9. 2D RDFs for phosphate−phosphate pairs with PSM for (A) POPC−PSM2:1, (B) POPC−PSM1:1, and (C) POPC−PSM1:2.

H bond

composition

POPE−PSM1:1 0.89 0.01 0.11 0.14

between POPC and PSM in POPC/PSM/cholesterol remains almost same as in POPC/PSM systems. POPE can also form intermolecular hydrogen bonds within itself, frequencies of which increase with increase in concentration of POPE. The average frequency of intermolecular hydrogen bond between POPE molecules in pure POPE system is 0.622 at 40 °C,61 where as it is 0.311 in system with 50% PSM in our study. On addition of PSM molecules, the hydrogen bond frequency 5203

DOI: 10.1021/acs.jpcb.7b00359 J. Phys. Chem. B 2017, 121, 5197−5208

Article

The Journal of Physical Chemistry B S10). For systems with cholesterol, the RDFs between cholesterol and the PSM (Figure 10) show two very sharp

Figure 11. Clustering of lipids for (A) POPC−PSM−cholesterol1:1:1 and (B) POPE−PSM−cholesterol1:1:1. Cholesterol: yellow; POPC/ POPE: red; PSM: blue.

motions, and τ3 represents longer motions (lipid axial), such as wobble and/or rotational motion. Table 6 contains the average time constants for the axial relaxation for C2−H and C2−C3 vectors for each system. Example correlation time plots of C2− C3 vector of POPC for the binary systems are shown in Figure S13. This large time constant is assumed to represent axial motion (lipid rotation and wobble motion).62 The other time constants are shown in Tables S4 and S5 and represent faster relaxing internal motions. An earlier simulation study with only PSM system,42 with same force field, showed τ1 = 0.559, τ2 = 20.41, and τ3 = 93.62 ns for C2−H. From our study, we see that internal motions represented by τ1 are slowed in the presence of POPC or POPE. The time constant for axial motion is 119.8 ± 25.4 ns in the presence of 67% PSM and 33% POPC and 63.65 ± 18.5 ns in the presence of 33% POPE. We see that with an increase in concentration of PSM the relaxation time τ3 increases, which is slowest in systems with cholesterol. Moreover, POPE at the low concentration results in a faster PSM axial relaxation compared to pure PSM bilayers due to altering the hydrogen bonding (or clustering) network. The slowest relaxation time is longer in systems with POPE than with POPC due to the hydrogen bonding of POPE with lipids. The internal motions are faster in the presence of more sphingomyelin. In the presence of cholesterol, the internal motions are slower in the system with POPE than with POPC.

Figure 10. 2D RDFs between cholesterol oxygen and carbonyl oxygen of POPC/POPE (blue) and carbonyl oxygen of PSM (red) for (A) POPE−PSM−Chol1:1:1and (B) POPC−PSM-Chol1:1:1.

peaks, i.e., first peak having g(r) value around 2.4 (Figure 10). A similar peak is observed between cholesterol and PC and PE lipids but at a g(r) that is less than 1, suggesting minor hydrogen bonding. To illustrate this, the hydrogen bond between the cholesterol oxygen and carbonyl oxygen of PSM in shown in Figure 8C. Overall, cholesterol prefers to cluster near the PSM lipid. Lipid clustering estimates the local packing of head groups in the bilayer because of hydrogen bonding and other lipid−lipid interactions. We see that average number of clusters do not change much with concentration (Table S3), but in the bilayers with POPE the number of clusters increase with increase in concentration of PSM. Number of clusters increase in the presence of cholesterol. Also in the systems with POPE, the number of clusters is slightly higher than the POPC systems. From the results, we see that head groups PC and PE dominate in the clusters more than PSM. In the systems with cholesterol, the number of cholesterol molecules is higher than the other two lipids. For POPE−PSM−cholesterol1:1:1, the snapshot (Figure 11B) shows that web-like connection of lipids (colored and those in gray) in the bilayer, suggesting an interaction network formed due to the intermolecular hydrogen bonds of POPE with PSM and cholesterol. Whereas in the system with POPC, there is clumped arrangement of lipids that lacks the web-like connection of lipids, as POPC forms less hydrogen bond with PSM. The plots for the binary systems are shown in Figures S11 and S12. Rotational or Axial Motions. Rotational relaxation motions provide details about lipid environment influencing lipid motions. From MD simulations, reorientational correlation functions are calculated for the C2−C3 (cross-chain) and C2−H vectors separately for POPC/POPE and PSM. The correlation functions (eq 2) are fitted to a relaxation model containing axial rotation or wobble and fast internal motions (eq 3).62 The time constants τ1 and τ2 represents the internal



DISCUSSION AND CONCLUSION In this work, we have studied the biophysical effects of sphingomyelin and cholesterols in the lipid bilayer containing PC or PE using all-atom molecular dynamics. Another objective was to develop simple eukaryotic membrane models containing sphingomyelin and cholesterol using accurate force fields. The simulations studies were performed using the latest all atom CHARMM force field which has been successfully tested on various lipid types.41,43,74 The presence of sphingomyelin and cholesterol introduces compactness in the bilayer, resulting in more ordered phase.73 The result is lowering of surface area. We have studied the effects at different concentrations of sphingomyelin in the bilayer. The SA/lipid provides information about the equilibrium of the simulations and provides information about the lipid packing and compactness.59,67 Our study shows that lipids tend to occupy less space in bilayer when the concentration of PSM increases. It is evident from the SA/lipid calculation where we find that the average surface area is higher in POPC−PSM and POPE−PSM bilayers for systems with 33% PSM concentration when compared to the systems with 66% PSM concentration (Table 1). The presence of cholesterol results in a more 5204

DOI: 10.1021/acs.jpcb.7b00359 J. Phys. Chem. B 2017, 121, 5197−5208

Article

The Journal of Physical Chemistry B

Table 6. Average Time Constants for the Axial Motion (in ns) of Three-Exponential Fitting to Reorientational Correlation Function with Standard Error of Mean for the C2−H and C2−C3 Vectors system POPC−PSM2:1 POPC−PSM1:1 POPC−PSM1:2 POPE−PSM2:1 POPE−PSM1:1 POPE−PSM1:2 POPC−PSM−cholesterol1:1:1 POPE−PSM−cholesterol1:1:1

PC/PE (C2−H) C2(t) τ3 59.0 74.8 56.5 51.5 37.3 126.8 138.5 169.6

± ± ± ± ± ± ± ±

PSM (C2−H) C2(t) τ3

17.0 22.1 11.4 18.5 8.2 87.2 44.3 23.6

48.8 86.2 119.8 44.6 35.73 63.7 94.3 128.1

± ± ± ± ± ± ± ±

10.6 20.6 25.4 0.3 3.2 18.5 18.6 22.3

PC/PE (C2−C3) C2(t) τ3 67.3 12.5 67.0 35.5 113.4 169.9 153.3 164.5

± ± ± ± ± ± ± ±

11.4 2.9 33.2 20.46 65.5 5.3 88.5 95.1

PSM (C2−C3) C2(t) τ3 46.7 83.8 130.1 36.6 31.8 95.4 119.9 127.9

± ± ± ± ± ± ± ±

7.6 26.2 27.7 2.5 4.5 25.5 44.4 19.4

POPC−PSM−Cholesterol1:1:1 system showed bilayer thickness of 44.0 ± 0.05 Å. The bilayer thickness for the pure POPC bilayer calculated by the same group was 36.9 Å. The bilayer thickness for pure POPC system calculated from experimental study was 37.9 Å,76 and that with CHARMM was 37.3 Å.61 Since CHARMM based bilayer thickness was closer to the experimental study, it is anticipated that for the system POPC− PSM−Cholesterol1:1:1 also the CHARMM-based bilayer thickness will be closer to the experimental value than the one based on Berger force field (united atom).18 Hydrogen bond analysis shows that the number of hydrogen bonds within PSM molecules and within POPE molecules increase with increase in respective concentrations. The detailed study of intra- and intermolecular hydrogen bonds within PSM molecules shows the frequency of hydrogen bond between O−H:::O−P is lower when compared to the pure PSM bilayers. The probability of the intermolecular hydrogen bond between amide nitrogen and phosphate oxygen increases, whereas that between amide nitrogen and carbonyl oxygen decreases from the pure PSM systems.42 Intermolecular hydrogen bond calculations between phospholipids and cholesterol in ternary systems show that sterols form more hydrogen bonds with sphingomyelin than phospholipids. This in agreement with the earlier study77 that sterols have a higher affinity for sphingomyelin than for phosphatidylcholine bilayers. In this study, we have built lipid bilayer models which could be used as simple models for studying the outer leaflet of eukaryotic plasma membranes rich in cholesterol and sphingomyelin. Lipid bilayers containing phosphatidylcholines, phosphatidylethanolamines, sphingomyelin, and cholesterol constitute model membranes for studying phase separation and domain formation in membranes.8 Composition and concentrations of lipid components in bilayers influence their activities and in turn influence the activities of associated proteins. Sphingomyelin-rich bilayer models may act as model membranes for studying proteins associated with myelin sheath, brain tissue, ocular lens, and erythrocytes. Sphingomyelin- and cholesterol-rich membranes are involved in many cellular processes such as membrane sorting, trafficking, signal transduction, cell polarization, and apoptosis.78 Accurate representation of different bilayers models will require incorporation of various lipids in different concentrations with proper force fields. The importance of sphingomyelin and cholesterol in development of more ordered lipid phase and lipid rafts was validated with overall decrease in SA/lipid and increase in order parameters in the model systems. In the binary mixtures, more sphingomyelin resulted in less SA/lipid and increase in order parameters and vice versa. The properties of the bilayer layer models built in this study agree with experimental studies with better representation using advanced

compact bilayer, as we see in our study, a large decrease in SA/ lipid from the systems without cholesterol (Table 1). Previous studies have been performed on this ternary system (POPC− PSM−Chol) with different concentrations of POPC33,73 using united atom force field and also with CHARMM force field.45 All these studies prove that there is a huge decrease in SA/lipid in the presence of high concentration of cholesterol and thus higher concentration of cholesterol results in more ordered phase. An earlier study (using Berger lipid force field and GROMOS parameters for headgroup)73 reported SA/lipid of 41 ± 1 Å for this particular ternary system with ratio POPC− PSM−Chol1:1:1. In our present study, for system POPC− PSM−Chol1:1:1, the SA/lipid was 44.4 ± 0.2 Å2. This difference may be due to different force field used in the study. The earlier study found the SA/lipid for pure POPC system to be 68.3 ± 0.01 Å using united atom force field. The SA/lipid for pure POPC (64.7 Å2) using CHARMM was in accordance to the experimentally determined values.61 A recent simulation study with the same CHARMM force field on a ternary mixture with the composition POPC−PSM−Chol0.08:0.61:0.31 (depicting liquid ordered phase) showed an average surface area of 40.5 ± 0.2 and 54.8 ± 0.5 Å for the composition POPC−PSM− Chol0.69:0.23:0.08 (depicting liquid disordered phase)45 which further validates that with increase in ordered phase of bilayer the SA/lipid decreases more. Deuterium order parameters provide information about the structural orientation and flexibility of lipids in bilayers.58 Order parameters calculated in our study are in excellent agreement with the experimental studies.39 The increase in order parameter with increase in PSM concentration corroborated well with the experimental studies.39,75 The increase in SCD values also supported the decrease in SA/lipid with increase in concentration of PSM. For systems, with cholesterol, the order parameters for the sn-1 chains of POPC and POPE were much higher than those of pure systems. The trend of order parameters along the carbon chain correlated better with pure systems when compared to those obtained from united atom based studies.33 Deuterium order parameters for sn-1 chain of POPE in POPE−PSM−Cholesterol1:1:1 correlated well with the experimental values (Figure S7). The order parameter for carbon atom 9′ in PSM acyl chain well correlated with recent NMR-based calculation71 for the ternary systems. The main discrepancy is a slightly higher SCD for lipids with binary mixture of POPC and PSM with MD simulations using the C36 lipid force field. This suggest that the surface area per lipid for these mixtures is slightly too low compared to experiments. The mass density profiles showed that the bilayer thickness increases with increase in concentration of PSM (Table 2) and the addition of cholesterol further increases the bilayer thickness. The earlier studies with united atom force field for 5205

DOI: 10.1021/acs.jpcb.7b00359 J. Phys. Chem. B 2017, 121, 5197−5208

Article

The Journal of Physical Chemistry B force field for lipids. Overall, this study justifies the change in biophysical properties of bilayers containing sphingomyelin and cholesterol along with phosphatidylcholines/phosphatidylethanolamines, and thus the bilayers can be used as models for studying various eukaryotic membranes.



(11) Lingwood, D.; Simons, K. Lipid rafts as a membrane-organizing principle. Science 2010, 327 (5961), 46−50. (12) Hjort Ipsen, J.; Karlstrom, G.; Mourtisen, O. G.; Wennerstrom, H.; Zuckermann, M. J. Phase-equilibria in the phosphatidylcholinecholesterol system. Biochim. Biophys. Acta, Biomembr. 1987, 905 (1), 162−172. (13) Leftin, A.; Job, C.; Beyer, K.; Brown, M. F. Solid-state 13C NMR reveals annealing of raft-like membranes containing cholesterol by the intrinsically disordered protein α-synuclein. J. Mol. Biol. 2013, 425 (16), 2973−2987. (14) Simons, K.; Ehehalt, R. Cholesterol, lipid rafts, and disease. J. Clin. Invest. 2002, 110 (5), 597−603. (15) Dickson, C. J.; Rosso, L.; Betz, R. M.; Walker, R. C.; Gould, I. R. GAFFlipid: A general amber force field for the accurate molecular dynamics simulation of phospholipid. Soft Matter 2012, 8 (37), 9617− 9627. (16) Skjevik, A. A.; Madej, B. D.; Walker, R. C.; Teigen, K. LIPID11: A modular framework for lipid simulations using amber. J. Phys. Chem. B 2012, 116 (36), 11124−11136. (17) Egberts, E.; Marrink, S. J.; Berendsen, H. J. C. Moleculardynamics simulation of a phospholipid membrane. Eur. Biophys. J. 1994, 22 (6), 423−436. (18) Berger, O.; Edholm, O.; Jahnig, F. Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature. Biophys. J. 1997, 72 (5), 2002−2013. (19) Holtje, M.; Forster, T.; Brandt, B.; Engels, T.; von Rybinski, W.; Holtje, H. D. Molecular dynamics simulations of stratum corneum lipid models: fatty acids and cholesterol. Biochim. Biophys. Acta, Biomembr. 2001, 1511 (1), 156−167. (20) Niemela, P.; Hyvonen, M. T.; Vattulainen, I. Structure and dynamics of sphingomyelin bilayer: Insight gained through systematic comparison to phosphatidylcholine. Biophys. J. 2004, 87 (5), 2976− 2989. (21) Mukhopadhyay, P.; Vogel, H. J.; Tieleman, D. P. Distribution of pentachlorophenol in phospholipid bilayers: A molecular dynamics study. Biophys. J. 2004, 86 (1), 337−345. (22) Mukhopadhyay, P.; Monticelli, L.; Tieleman, D. P. Molecular dynamics simulation of a palmitoyl-oleoyl phosphatidylserine bilayer with Na+ counterions and NaCl. Biophys. J. 2004, 86 (3), 1601−1609. (23) MacDermaid, C. M.; Kashyap, H. K.; DeVane, R. H.; Shinoda, W.; Klauda, J. B.; Klein, M. L.; Fiorin, G. Molecular dynamics simulations of cholesterol-rich membranes using a coarse-grained force field for cyclic alkanes. J. Chem. Phys. 2015, 143, 243144. (24) Monje-Galvan, V.; Klauda, J. B. Modeling yeast organelle membranes and how lipid diversity influences bilayer properties. Biochemistry 2015, 54 (45), 6852−6861. (25) Khakbaz, P.; Klauda, J. B. Probing the importance of lipid diversity in cell membranes via molecular simulation. Chem. Phys. Lipids 2015, 192, 12−22. (26) Pandit, K. R.; Klauda, J. B. Membrane models of e-coli containing cyclic moieties in the aliphatic lipid chain. Biochim. Biophys. Acta, Biomembr. 2012, 1818 (5), 1205−1210. (27) O’Connor, J. W.; Klauda, J. B. Lipid membranes with a majority of cholesterol: Applications to the ocular lens and aquaporin 0. J. Phys. Chem. B 2011, 115 (20), 6455−6464. (28) Ingolfsson, H. I.; Melo, M. N.; van Eerden, F. J.; Arnarez, C.; Lopez, C. A.; Wassenaar, T. A.; Periole, X.; de Vries, A. H.; Tieleman, D. P.; Marrink, S. J. Lipid organization of the plasma membrane. J. Am. Chem. Soc. 2014, 136 (41), 14554−14559. (29) Pan, J. J.; Cheng, X. L.; Monticelli, L.; Heberle, F. A.; Kucerka, N.; Tieleman, D. P.; Katsaras, J. The molecular structure of a phosphatidylserine bilayer determined by scattering and molecular dynamics simulations. Soft Matter 2014, 10 (21), 3716−3725. (30) Bennett, W. F. D.; Tieleman, D. P. Molecular simulation of rapid translocation of cholesterol, diacylglycerol, and ceramide in model raft and nonraft membranes. J. Lipid Res. 2012, 53 (3), 421− 429.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b00359. Figures for SA/lipid for all systems, order parameters comparing simulation with experiment when available and/or comparing with pure lipid system, component mass density profiles for all the systems, 2D RDFs for phosphate−phosphate pairs, clustering plots for binary systems and time correlation function, C(t), vs time for C2−C3 vector; tables with details of the bilayer systems, peak values of RDFs, lipid clustering results, average time constants for the C2−H vector (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected], Ph (301) 405-1320 (J.B.K.). ORCID

Jeffery B. Klauda: 0000-0001-8725-1870 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This research is supported by NSF Grant MCB-1149187. Simulations were performed on High Performance Computing (HPC) resource (Deepthought1 and Deepthought2) of Division of Information Technology at University of Maryland, College Park and HPC of Maryland Advanced Research Computing Center (MARCC) managed by John Hopkins University and University of Maryland. We thank Xiaohong Zhuang at UMD for her help in calculations of component SA and lipid clustering.



REFERENCES

(1) van Meer, G.; Voelker, D. R.; Feigenson, G. W. Membrane lipids: Where they are and how they behave. Nat. Rev. Mol. Cell Biol. 2008, 9 (2), 112−124. (2) Alberts, B. Molecular Biology of the Cell, 5th ed.; Garland Science: New York, 2008. (3) Opdenkamp, J. A. F. Lipid asymmetry in membranes. Annu. Rev. Biochem. 1979, 48, 47−71. (4) Jacobson, K.; Sheets, E. D.; Simson, R. Revisiting the fluid mosaic model of membranes. Science 1995, 268 (5216), 1441−1442. (5) Brown, D. A.; London, E. Structure and function of sphingolipidand cholesterol-rich membrane rafts. J. Biol. Chem. 2000, 275 (23), 17221−17224. (6) Heberle, F. A.; Feigenson, G. W. Phase separation in lipid membranes. Cold Spring Harbor Perspect. Biol. 2011, 3, a004630. (7) de Almeida, R. F. M.; Fedorov, A.; Prieto, M. Sphingomyelin/ phosphatidylcholine/cholesterol phase diagram: Boundaries and composition of lipid rafts. Biophys. J. 2003, 85 (4), 2406−2416. (8) Edidin, M. The state of lipid rafts: From model membranes to cells. Annu. Rev. Biophys. Biomol. Struct. 2003, 32, 257−283. (9) Simons, K.; Vaz, W. L. C. Model systems, lipid rafts, and cell membranes. Annu. Rev. Biophys. Biomol. Struct. 2004, 33, 269−295. (10) Silvius, J. R. Role of cholesterol in lipid raft formation:Lessons from lipid model systems. Biochim. Biophys. Acta, Biomembr. 2003, 1610 (2), 174−183. 5206

DOI: 10.1021/acs.jpcb.7b00359 J. Phys. Chem. B 2017, 121, 5197−5208

Article

The Journal of Physical Chemistry B (31) Zhang, Z.; Bhide, S. Y.; Berkowitz, M. L. Molecular dynamics simulations of bilayers containing mixtures of sphingomyelin with cholesterol and phosphatidylcholine with cholesterol. J. Phys. Chem. B 2007, 111 (44), 12888−12897. (32) Khelashvili, G. A.; Scott, H. L. Combined monte carlo and molecular dynamics simulation of hydrated 18:0 sphingomyelincholesterol lipid bilayers. J. Chem. Phys. 2004, 120 (20), 9841−9847. (33) Aittoniemi, J.; Niemela, P. S.; Hyvonen, M. T.; Karttunen, M.; Vattulainen, I. Insight into the putative specific interactions between cholesterol, sphingomyelin, and palmitoyl-oleoyl phosphatidylcholine. Biophys. J. 2007, 92 (4), 1125−1137. (34) Waheed, Q.; Tjornhammar, R.; Edholm, O. Phase transitions in coarse-grained lipid bilayers containing cholesterol by molecular dynamics simulations. Biophys. J. 2012, 103 (10), 2125−2133. (35) Hofsass, C.; Lindahl, E.; Edholm, O. Molecular dynamics simulations of phospholipid bilayers with cholesterol. Biophys. J. 2003, 84 (4), 2192−2206. (36) Rog, T.; Pasenkiewicz-Gierula, M.; Vattulainen, I.; Karttunen, M. Ordering effects of cholesterol and its analogues. Biochim. Biophys. Acta, Biomembr. 2009, 1788 (1), 97−121. (37) Niemela, P. S.; Hyvonen, M. T.; Vattulainen, I. Atom-scale molecular interactions in lipid raft mixtures. Biochim. Biophys. Acta, Biomembr. 2009, 1788 (1), 122−135. (38) Chiu, S. W.; Vasudevan, S.; Jakobsson, E.; Mashl, R. J.; Scott, H. L. Structure of sphingomyelin bilayers: A simulation study. Biophys. J. 2003, 85 (6), 3624−3635. (39) Mehnert, T.; Jacob, K.; Bittman, R.; Beyer, K. Structure and lipid interaction of N-palmitoylsphingomyelin in bilayer membranes as revealed by H-2-NMR spectroscopy. Biophys. J. 2006, 90 (3), 939− 946. (40) Bartels, T.; Lankalapalli, R. S.; Bittman, R.; Beyer, K.; Brown, M. F. Raftlike mixtures of sphingomyelin and cholesterol investigated by solid-state H-2 NMR spectroscopy. J. Am. Chem. Soc. 2008, 130 (44), 14521−14532. (41) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D.; Pastor, R. W. Update of the CHARMM all-atom additive force field for lipids: Validation on six lipid types. J. Phys. Chem. B 2010, 114 (23), 7830−7843. (42) Venable, R. M.; Sodt, A. J.; Rogaski, B.; Rui, H.; Hatcher, E.; MacKerell, A. D.; Pastor, R. W.; Klauda, J. B. CHARMM all-atom additive force field for sphingomyelin: Elucidation of hydrogen bonding and of positive curvature. Biophys. J. 2014, 107 (1), 134−145. (43) Lim, J. B.; Rogaski, B.; Klauda, J. B. Update of the cholesterol force field parameters in CHARMM. J. Phys. Chem. B 2012, 116 (1), 203−210. (44) Pluhackova, K.; Kirsch, S. A.; Han, J.; Sun, L. P.; Jiang, Z. Y.; Unruh, T.; Bockmann, R. A. A critical comparison of biomembrane force fields: Structure and dynamics of model DMPC, POPC, and POPE bilayers. J. Phys. Chem. B 2016, 120 (16), 3888−3903. (45) Sodt, A. J.; Pastor, R. W.; Lyman, E. Hexagonal substructure and hydrogen bonding in liquid-ordered phases containing palmitoyl sphingomyelin. Biophys. J. 2015, 109 (5), 948−955. (46) Jo, S.; Lim, J. B.; Klauda, J. B.; Im, W. CHARMM-GUI membrane builder for mixed bilayers and its application to yeast membranes. Biophys. J. 2009, 97 (1), 50−58. (47) Jo, S.; Kim, T.; Iyer, V. G.; Im, W. Software news and updates CHARMM-GUI: A web-based graphical user interface for CHARMM. J. Comput. Chem. 2008, 29 (11), 1859−1865. (48) Jo, S.; Kim, T.; Im, W. Automated builder and database of protein/membrane complexes for molecular dynamics simulations. PLoS One 2007, 2 (9), e880. (49) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79 (2), 926−935. (50) Durell, S. R.; Brooks, B. R.; Ben-Naim, A. Solvent-induced forces between two hydrophilic groups. J. Phys. Chem. 1994, 98 (8), 2198−2202.

(51) Lee, J.; Cheng, X.; Jo, S.; MacKerell, A. D.; Klauda, J. B.; Im, W. CHARMM-GUI input generator for NAMD, Gromacs, Amber, Openmm, and CHARMM/OpenMM simulations using the CHARMM36 additive force field. Biophys. J. 2016, 110 (3), 641a. (52) Vandrunen, R.; Vanderspoel, D.; Berendsen, H. J. C. Gromacs a software package and a parallel computer for molecular-dynamics. Abstr. Pap. Am. Chem. Soc. 1995, 209, 49-COMP. (53) Hoover, W. G. Canonical Dynamics - Equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31 (3), 1695− 1697. (54) Parrinello, M.; Rahman, A. Polymorphic transitions in singlecrystals - a new molecular-dynamics method. J. Appl. Phys. 1981, 52 (12), 7182−7190. (55) Darden, T.; York, D.; Pedersen, L. Particle mesh ewald - an n.log(n) method for ewald sums in large systems. J. Chem. Phys. 1993, 98 (12), 10089−10092. (56) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18 (12), 1463−1472. (57) Barber, C. B.; Dobkin, D. P.; Huhdanpaa, H. The Quickhull algorithm for convex hulls. ACM Trans. Math. Software 1996, 22 (4), 469−483. (58) Leftin, A.; Brown, M. F. An NMR database for simulations of membrane dynamics. Biochim. Biophys. Acta, Biomembr. 2011, 1808 (3), 818−839. (59) Petrache, H. I.; Dodd, S. W.; Brown, M. F. Area per lipid and acyl length distributions in fluid phosphatidylcholines determined by H-2 NMR spectroscopy. Biophys. J. 2000, 79 (6), 3172−3192. (60) Schindler, H.; Seelig, J. Deuterium order parameters in relation to thermodynamic properties of a phospholipid bilayer - statistical mechanical interpretation. Biochemistry 1975, 14 (11), 2283−2287. (61) Zhuang, X.; Dávila-Contreras, E. M.; Beaven, A. H.; Im, W.; Klauda, J. B. An extensive simulation study of lipid bilayer properties with different head groups, acyl chain lengths, and chain saturations. Biochim. Biophys. Acta, Biomembr. 2016, 1858 (12), 3093−3104. (62) Klauda, J. B.; Roberts, M. F.; Redfield, A. G.; Brooks, B. R.; Pastor, R. W. Rotation of lipids in membranes: Molecular dynamics simulation, P-31 spin-lattice relaxation, and rigid-body dynamics. Biophys. J. 2008, 94 (8), 3074−3083. (63) Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; Vanderplas, J.; Passos, A.; Cournapeau, D.; Brucher, M.; Perrot, M.; Duchesnay, E. Scikit-learn: Machine learning in python. J. Mach. Learn. Res. 2011, 12, 2825−2830. (64) Sander, J.; Ester, M.; Kriegel, H. P.; Xu, X. W. Density-based clustering in spatial databases: The algorithm GDBSCAN and its applications. Data Min. Knowl. Discovery 1998, 2 (2), 169−194. (65) Ester, M.; Kriegel, H. P.; Sander, J.; Xu, X. In A density-based algorithm for discovering clusters in large spatial databases with noise, Proceedings of 2nd International Conference on knowledge discovery and data mining, AAAI Press: 1996; pp 226−231. (66) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14 (1), 33−38. (67) Kinnun, J. J.; Mallikarjunaiah, K. J.; Petrache, H. I.; Brown, M. F. Elastic deformation and area per lipid of membranes: Atomistic view from solid-state deuterium NMR spectroscopy. Biochim. Biophys. Acta, Biomembr. 2015, 1848 (1), 246−259. (68) de Meyer, F.; Smit, B. Effect of cholesterol on the structure of a phospholipid bilayer. Proc. Natl. Acad. Sci. U. S. A. 2009, 106 (10), 3654−3658. (69) Bunge, A.; Muller, P.; Stockl, M.; Herrmann, A.; Huster, D. Characterization of the ternary mixture of sphingomyelin, POPC, and cholesterol: Support for an inhomogeneous lipid distribution at high temperatures. Biophys. J. 2008, 94 (7), 2680−2690. (70) Kucerka, N.; van Oosten, B.; Pan, J. J.; Heberle, F. A.; Harroun, T. A.; Katsaras, J. Molecular structures of fluid phosphatidylethanolamine bilayers obtained from simulation-to-experiment comparisons and experimental scattering density profiles. J. Phys. Chem. B 2015, 119 (5), 1947−1956. 5207

DOI: 10.1021/acs.jpcb.7b00359 J. Phys. Chem. B 2017, 121, 5197−5208

Article

The Journal of Physical Chemistry B (71) Engberg, O.; Yasuda, T.; Hautala, V.; Matsumori, N.; Nyholm, T. K. M.; Murata, M.; Slotte, J. P. Lipid interactions and organization in complex bilayer membranes. Biophys. J. 2016, 110 (7), 1563−1573. (72) Zhuang, X. H.; Makover, J. R.; Im, W.; Klauda, J. B. A systematic molecular dynamics simulation study of temperature dependent bilayer structural properties. Biochim. Biophys. Acta, Biomembr. 2014, 1838 (10), 2520−2529. (73) Niemela, P. S.; Ollila, S.; Hyvonen, M. T.; Karttunen, M.; Vattulainen, I. Assessing the nature of lipid raft membranes. PLoS Comput. Biol. 2007, 3 (2), e34. (74) Lim, J. B.; Klauda, J. B. Lipid chain branching at the iso- and anteiso-positions in complex chlamydia membranes: A molecular dynamics study. Biochim. Biophys. Acta, Biomembr. 2011, 1808 (1), 323−331. (75) Soni, S. P.; LoCascio, D. S.; Liu, Y. D.; Williams, J. A.; Bittman, R.; Stillwell, W.; Wassall, S. R. Docosahexaenoic acid enhances segregation of lipids between raft and nonraft domains: H-2-NMR study. Biophys. J. 2008, 95 (1), 203−214. (76) Kucerka, N.; Tristram-Nagle, S.; Nagle, J. F. Structure of fully hydrated fluid phase lipid bilayers with monounsaturated chains. J. Membr. Biol. 2006, 208 (3), 193−202. (77) Lonnfors, M.; Doux, J. P. F.; Killian, J. A.; Nyholm, T. K. M.; Slotte, P. Sterols have higher affinity for sphingomyelin than for phosphatidylcholine bilayers even at equal acyl-chain order. Biophys. J. 2011, 100 (11), 2633−2641. (78) Giocondi, M. C.; Boichot, S.; Plenat, T.; Le Grimellec, C. Structural diversity of sphingomyelin microdomains. Ultramicroscopy 2004, 100 (3−4), 135−143.

5208

DOI: 10.1021/acs.jpcb.7b00359 J. Phys. Chem. B 2017, 121, 5197−5208