Influence of Cholesterol on Phospholipid Bilayer Structure and

Oct 22, 2016 - (22) NAMD software package(23) was used to run 130-ns trajectories of each model to ensure proper equilibration. ...... Pan , J.; Marqu...
0 downloads 6 Views 8MB Size
Download PDF
Article pubs.acs.org/JPCB

Influence of Cholesterol on Phospholipid Bilayer Structure and Dynamics Christopher T. Boughter,† Viviana Monje-Galvan,† Wonpil Im,‡ and Jeffery B. Klauda*,†,§ †

Department of Chemical and Biomolecular Engineering, University of Maryland, College Park, Maryland 20742, United States; Department of Biological Sciences and Bioengineering Program, Lehigh University, Bethlehem, Pennsylvania 18015, United States § Biophysics Program, University of Maryland, College Park, Maryland 20742, United States ‡

S Supporting Information *

ABSTRACT: In this study, the influence of cholesterol on lipid bilayers is investigated by changing phospholipid headgroup, cholesterol concentration, chain saturation, and temperature. Molecular dynamics (MD) simulations were used to characterize bilayers containing phosphatidylcholine (PC) head groups with either fully saturated dimyristoyl (DM) or monounsaturated dioleoyl (DO) acyl chains and cholesterol concentrations ranging from 5 to 50%. To further explore the effects of cholesterol on bilayers with different head groups, we also performed MD simulations of bilayer systems having 15% cholesterol with phosphatidic acid (PA), phosphatidylethanolamine (PE), phosphatidylglycerol (PG), phosphatidylinositol (PI), and phosphatidylserine (PS), each having DM chains and at a temperature above the solid gel phase transition. Additionally, bilayers of DMPA, DMPE, and DMPS with 15% cholesterol were simulated at temperatures below the solid gel phase transition temperatures. Compared to membranes without cholesterol, cholesterol in the model bilayers increases chain order in bilayers with the highest order in the liquid ordered and solid gel phases. Head group properties and acyl chain saturation are also found to critically impact bilayer dynamics, largely through the formation of hydrogen bonds between membrane components. These results provide a better understanding of the basic characteristics on structure and dynamics of cholesterol-containing membranes by revealing molecular details of interactions between cholesterol and phospholipids as well as add to the library of simulation data necessary for the MD community to accurately represent relevant models of atomic-scale systems.

■

INTRODUCTION Surrounding every living cell is a membrane comprised of an assortment of proteins and lipids arranged in a bilayer. These bilayers are dynamic in nature with lipid molecules diffusing within a leaflet and occasionally flipping between leaflets.1 The phases of such dynamic systems have been intensely studied, yet still are not well understood.2 The different phases of a bilayer are important for the proper functionality of membrane proteins and the cell as a whole.3 Membrane proteins can be either anchored to the membrane surface, superficially attached to the membrane, or traverse the membrane completely.4 A bilayer that is too fluid or too rigid would alter interactions between its components as well as with proteins, possibly leading to reduced protein functionality. Among the proteins whose functions are altered by membrane fluidity are those involved in signal transduction, as their proper function is dependent on mobility within the membrane.5 Therefore, the ability to alter membrane fluidity is vital to cells to maintain homeostasis in environments with frequent temperature fluctuations,6 and these vital changes are mediated by altering the concentration of lipids in bilayers.7 The phospholipids used in this study are all critical to important functions within the cell. Phosphatidic acid (PA) is an important intermediary for the synthesis of glycerophos© 2016 American Chemical Society

pholipids and acts as a determining factor in the curvature of a membrane.8,9 Phosphatidylcholine (PC), phosphatidylethanolamine (PE), and phosphatidylinositol (PI) all have the potential to be hydrolyzed and used to propagate signals through their respective pathways.10 Phosphatidylglycerol (PG) is one of the most prolific anionic lipids in bacterial cells and is critical to cell viability.11 Lastly, phosphatidylserine (PS) plays a crucial role in signaling phagocytes to remove dead or dying cells during apoptosis.12 While head groups play an important role in phospholipid interactions, saturated and monounsaturated acyl chains can alter bilayer dynamics in significant ways. Phospholipids with saturated acyl chains tend to condense a bilayer and result in higher melting temperatures than bilayers containing unsaturated chains, primarily due to the kink in the unsaturated carbon−carbon cis double bond.13 Cholesterol, one of many sterols present in cell membranes of living organisms, is the most common in animal cells.14 Cholesterol is comprised of steroidal rings with an additional covalently bonded hydroxyl group (Figure S1). Previous studies have shown that when inserted into a phospholipid bilayer, Received: August 24, 2016 Revised: October 19, 2016 Published: October 22, 2016 11761

DOI: 10.1021/acs.jpcb.6b08574 J. Phys. Chem. B 2016, 120, 11761−11772

Article

The Journal of Physical Chemistry B

Table 1. Summary of Simulated Models in This Studya

cholesterol acts as a normalizing agent, creating disorder in highly ordered bilayers and bringing order to those that are disordered.15 Over the years there has been an ever-increasing interest in the dynamics of these effects, particularly in the search for the formation of lipid rafts. Recent results suggest that these rafts develop from repulsive interactions between cholesterol and phospholipids which melt at relatively lower temperatures.3 Regardless of the validity of its often-debated biophysical role, cholesterol is a critical contributor to bilayer dynamics in all phases. Molecular dynamics (MD) simulations provide a microscope through which researchers can peer into the atomistic interactions in membranes, which are difficult to obtain from experiments. Here we report MD simulations to explore atomistic interactions between cholesterol and various phospholipids to gain insight into the physical mechanisms behind some of the larger scale interactions in the various bilayers that encapsulate cells and organelles in living systems. This study contributes to the overall effort to understand the structural and dynamical phase behavior of phospholipid bilayers with cholesterol.

a

lipid

% cholesterol

temperature (K)

# systems

DMPC DOPC DMPA DMPE DMPS DMPG DMPI

5, 15, 20, 40, 50 5, 15, 20, 30, 40, 50 15 15 15 15 15

303.15 303.15 303.15, 333.15 303.15, 333.15 303.15, 333.15 303.15 303.15

5 6 2 2 2 1 1

See Table S1 for a detailed description of each system.

dioleoyl (DO) 16-carbon monounsaturated chain. To further probe bilayer dynamics and structure, simulations were run for DM chains in bilayers having a 15% cholesterol concentration with different lipid head groups, including PA, PE, PG, PI, and PS. As a final test DMPA, DMPE, and DMPS bilayers were also simulated with 15% cholesterol content at 333.15 K, i.e., above their pure gel-phase transition temperatures. Throughout the article, the mention of this gel-phase transition temperature is meant to serve as a reference to how the pure phospholipid bilayer behaves. For bilayers with cholesterol, the solid gel phase is thermodynamically stable at the appropriate composition and temperature.15 Simulation Analysis. After equilibration of the systems, the trajectories were analyzed for direct comparison between simulations and available experimental data. We performed statistical analysis on the overall and component-wise surface area per lipid (SA), deuterium order parameters (SCD), lipid wobble (slow rotation) time scales, sterol tilt angles, radial distribution functions (RDFs), and electron density profiles (EDPs) for all models. Surface Area per Lipid (SA). Stability of the simulation was determined from the overall SA for each trajectory. The total SA was calculated from the XY plane of the membrane and divided by the number of lipids per leaflet. In addition to being a good indicator of simulation equilibration, the SA can also provide insight into lipid packing of the membrane as well as the impact that cholesterol has on bilayers. Individual component SA was computed using Qhull,33 a software that uses a triangulation technique to generate a polygon from the coordinates of three atoms in a lipid (C2-carbon connecting sn2 acyl chain, C31-first carbon on sn-1, C21-first carbon on sn2), or one atom in a sterol (O3-hydroxyl oxygen), and determines its volume and surface area in a given plane. Deuterium Order Parameters (SCD). NMR SCD values are a measure of the orientational mobility of the carbon−deuterium (C−D) bond and specifically measure the C−D orientation with respect to the bilayer normal. An SCD value of zero corresponds to either a completely disordered system or a perfectly ordered system at an angle of 54.7° with respect to the bilayer normal.34 These segmental parameters measure the relative order inside a bilayer. In simulations, the C−H vector is used to obtain SCD values that can be directly compared to NMR experimental data.7 SCD values were computed for each model bilayer system using eq 1,

■

METHODS MD simulations are necessary to study changes in a system over nanosecond to microsecond time scales in atomistic detail. Twenty-one membrane models were developed to study the impact of interactions between head groups, chain saturation, and cholesterol on bilayer structure and the dynamics of phase behavior. Each simulated model was built using CHARMMGUI Membrane Builder,16−20 contained 100 lipids evenly distributed in two leaflets with neutralizing ions, and was fully hydrated using the TIP3P model for water (see Table S1 for system setup).21 All models were equilibrated on CHARMM using the standard six-step CHARMM-GUI protocol for 225 ps and periodic boundary conditions.16,18 The first two steps were run using NVT (constant particle number, volume, and temperature) dynamics, and the remaining four using the NPT (constant particle number, pressure, and temperature) ensemble to slowly reduce various restraints on the system.22 NAMD software package23 was used to run 130-ns trajectories of each model to ensure proper equilibration. We ran three replicates of each model with the CHARMM36 lipid force field,24 which includes the most updated and accurate parameters for PI lipids25 and sterols26 needed for this study. The NPT ensemble was used during all production runs with a 2 fs time-step and the SHAKE algorithm to constrain the bonds having hydrogen atoms.27 The temperature was kept constant at 303.15 K using the Nosé−Hoover thermostat and Langevin dynamics, with the exception of three model membranes which were simulated at 333.15 K to observe temperature influence on bilayer structure and phase dynamics (Table 1).28 For constant pressure control at 1 bar we used the Langevin piston on CHARMM and the Nosé−Hoover Langevin piston on NAMD allowing the cell box size to change semi-isotropically.29,30 van der Waals interactions were computed using a Lennard-Jones force-switching function over 10 to 12 Å, while long-range electrostatic interactions used particle mesh Ewald.31,32 Overall, three replicates of twenty-one distinct systems were run for a total of 7.67 μs of simulation time. Initial simulations looked at bilayers with varying cholesterol concentrations ranging from 5 to 50%. PC head groups were simulated with a dimyristoyl (DM) 14-carbon saturated fatty acid chain or a

SCD =

3 2 1 cos β − 2 2

(1)

where β is the angle between a given C−H bond in the lipid tail and the bilayer normal. Each tail was defined by its position with respect to the glycerol group. By convention, the sn-2 chain is defined as the tail attached to the oxygen atom of the 11762

DOI: 10.1021/acs.jpcb.6b08574 J. Phys. Chem. B 2016, 120, 11761−11772

Article

The Journal of Physical Chemistry B

Figure 1. Snapshots of the DMPC (A,B) and DMPI (C,D) bilayers at 303.15 K with 15% cholesterol. DMPC is colored red, DMPI is yellow, cholesterol is black, and water is blue. Dotted black box indicates the primary simulation system, with everything outside the box corresponding to image atoms from the periodic boundary conditions.

of each cholesterol molecule over the equilibrated trajectory, which allows us to quantify cholesterol-cholesterol interactions on the bilayer plane. Hydrogen Bonding Characterization. The hbond functionality in CHARMM program39 was used to count the hydrogen bonds between lipid molecules, both sterol-lipid and lipid− lipid. This was done including image atoms near the periodic boundary to take into account the hydrogen bonds of the molecules at the border of the main simulation box. The options for the hbond command were set to 2.4 Å and 150° for the cutoff hydrogen-heavy atom distance and angle, respectively. The values reported in this study are averaged across the equilibrated trajectory of all three replicates for each system. Electron Density Profiles (EDPs). The structure of biological lipid membranes is difficult to obtain experimentally due to thermal fluctuations. Neutron and X-ray scattering experiments such as small angle neutron scattering (SANS) and small-angle X-ray scattering (SAXS) can be used to obtain an average of these fluctuations to give a general idea of a bilayer structure.40,41 EDPs from simulation data were processed using the SIMtoEXP software package following the procedure of Kučerka et al.42 to determine the bilayer head-to-head distance (DHH), the Luzzati thickness (overall bilayer thickness: DB), and its hydrocarbon region or core (2DC).42 DB and 2DC were estimated from 50% of the volume probability distributions of water and fatty acid tails, respectively. The EDPs were also compared to X-ray scattering experiments in Fourier space to determine simulation accuracy.

second carbon of the glycerol group; while the other tail is denoted as the sn-1 chain. The hydrogen atoms attached to the second carbon of the sn-2 chain give two different experimental signals, which were also computed from the simulation data.35−37 Lipid Tilt. To complete the molecular structural analysis of each system we computed the average sterol tilt angle as described by Lim et al.26 The tilt angle was defined as the angle between the bilayer normal and the vector defined between C17−C3, the carbons at the base of the sterol and the one attached to its hydroxyl group, respectively. Similarly, lipid tail tilt was evaluated for some systems defining the angle between the bilayer normal and the vector between the last and first carbons in the tail. The tilt angle was computed for each tail separately and then averaged. Lipid Relaxation Times. We characterized the lipid motions using reorientational correlation functions and focused on axial relaxation time constants, which are relevant to determine the environment inside a bilayer from simulation data. These relaxation times were computed using the second rank reorientational correlation function C2(t), eq 2, for specific atoms in the simulation data: ̂ ·μ(̂ t ))⟩ C2(t ) = ⟨P2(μ(0)

(2)

where P2 is the second Legendre polynomial. μ̂ is a heavy atomhydrogen vector selected to estimate the slow relaxation time (classified as wobble38) between the first carbon after the carbonyl of each fatty acid tail (C22 and C32) over at least 80 ns of equilibrated trajectory. Three-exponential fits were used to determine the time constants for each lipid in the models; eq 3 shows the general form of the fitting function.

■

RESULTS AND DISCUSSION All systems were analyzed to examine the effects of cholesterol concentration and lipid headgroup on bilayer structure. All reported values are the average over three replicas for each system and include their respective standard errors. The SA was calculated for each system as a measure of equilibrium and to gain insight into the lateral packing of lipids in the membrane. Acyl chain SCD were calculated to ascertain the ordering of lipids inside each bilayer and compare with available experimental data. Additionally, RDFs and hydrogen bonding frequencies were computed to determine local organization of lipids. EDPs were used to determine the position of lipid tails, phosphate, and headgroups with respect to the center of the bilayer, and estimate bilayer thicknesses. To look at molecular-

3

C2(t ) = ao +

∑ a i e −t / τ

i

i=1

(3)

The independent coefficient ao in eq 3 was estimated from the average value of the plateau reached by the correlation function over 15 to 45 ns (depending on the system). Radial Distribution Functions (RDFs). RDFs measure how the density of a given atom varies as a function of distance from a reference atom, and can be used to improve our understanding of the lateral ordering of bilayer components. RDFs were calculated between the cholesterol hydroxyl oxygen 11763

DOI: 10.1021/acs.jpcb.6b08574 J. Phys. Chem. B 2016, 120, 11761−11772

Article

The Journal of Physical Chemistry B

transition temperature.48 Varying lipid headgroup in systems with fully saturated tails (DM) resulted in similar SA values when simulated at temperatures above the corresponding gel transition temperature (Table 2) suggesting that headgroups do not play a large role in determining the surface area of the bilayer (in the presence of 15% cholesterol content), in agreement with experimental and computational results.49,50 However, there is a general trend in the overall and component SA that range from 48.9 (DMPI) to 52.4 Å2 (DMPG) in the following order: DMPI < DMPC < DMPS = DMPE < DMPA < DMPG. As expected, the decrease in temperature that takes PA, PE, and PS to the solid gel phase leads to a decrease in total area, as the acyl chains pack together more tightly and condense the bilayer. Additional MD simulations were run for DOPC (50% cholesterol) and DMPC (40% cholesterol) with membrane models having 100 lipids per leaflet to examine the effect of system size on the measured properties for liquid ordered (Lo) bilayers. There was no statistically significant difference in the SA of the small systems compared to the larger ones (see dopc50b and dmpc40b in Table S2), i.e., the values agree within 95% confidence interval (CI). Figure S2 shows the SCD values and cholesterol-cholesterol 2D RDFs plots comparing the data from the trajectories of the large and small systems. As expected, agreement in the SA also results in agreement with other structural membrane properties.51 In addition, lipid slow rotation time constants (wobble) showed no significant difference at 95% CI (see Table 3). Therefore, our use of 50lipid-per-leaflet systems is appropriate for the scope of this article.

level dynamical motions, slow relaxation of wobble motion were estimated with exponential fits to the reorientational correlation function. Sample end snapshot graphics were created in VMD43 and are shown in Figure 1. Lipid Lateral Packing. Figure 2A shows that there is a nearly linear decrease in the SA for DOPC as cholesterol

Table 3. Lipid Wobble (Slow Rotation) Times for All Lipids in This Study with ± Standard Errors

Figure 2. (A) Surface area per lipid (i.e., the total area divided by the total number of lipids in one leaflet) and (B) surface area per phospholipid (excluding cholesterol) of DMPC and DOPC bilayers at 303.15 K as a function of cholesterol concentration. The DMPC systems show the nonuniform condensing of the bilayer with increasing cholesterol concentration.

concentration increases, indicating that the system behaves as an ideal solution in agreement with experimental trends that show decrease in SA as cholesterol content increases.44−46 Conversely, the decrease is clearly nonlinear for DMPC, suggesting nonideal solution (i.e., nonrandom mixing). The same trend is observed for the component-wise SA in Figure 2B and Table S2. The lateral packing for DOPC is less than that of DMPC for all cholesterol concentrations. This is due in part to the kinked monounsaturated tails in DOPC. DMPA, DMPE, and DMPS bilayers have pure gel transition temperatures of 325, 323, and 308 K, respectively. 47 Simulations were carried at 333.15 K to show the liquiddisordered (Ld) behavior of these bilayers, and at 303.15 K to show the solid gel behavior. The simulation temperature was set to 303.15 K for DMPG and DMPI systems since that temperature is already above the corresponding pure gel

system

chol%

temp (K)

τ3 (ns)

DMPC5 DMPC15 DMPC20 DMPC40 DMPC40b DMPC50 DOPC5 DOPC15 DOPC20 DOPC30 DOPC40 DOPC50 DOPC50b DMPA DMPE DMPS DMPG DMPI DMPA2 DMPE2 DMPS2

0.05 0.15 0.2 0.4 0.4 0.5 0.05 0.15 0.2 0.3 0.4 0.5 0.4 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15

303.15

10.1 ± 1.9 12.7 ± 1.1 14.2 ± 1.1 24.4 ± 5.1 33.0 ± 4.6 22.0 ± 5.2 9.9 ± 4.0 9.9 ± 4.4 7.9 ± 0.6 10.7 ± 2.0 12.3 ± 1.6 13.8 ± 1.0 19.9 ± 4.6 11.9 ± 5.1 20.1 ± 1.1 23.9 ± 5.4 7.1 ± 0.1 29.7 ± 8.6 3.9 ± 1.3 5.1 ± 2.3 7.4 ± 1.0

333.15

Table 2. Average SA per Lipid [Å2] of the Systems with Varying Headgroup and 15% Cholesterol Content with ± Standard Errors system

DMPA

DMPE

DMPG

DMPI

DMPS

DMPC

liquid disordered solid gel

51.8 ± 0.1 42.6 ± 0.1

50.1 ± 0.1 42.5 ± 0.1

52.4 ± 0.1

48.9 ± 0.1

50.9 ± 0.1 43.8 ± 0.1

49.5 ± 0.1

11764

DOI: 10.1021/acs.jpcb.6b08574 J. Phys. Chem. B 2016, 120, 11761−11772

Article

The Journal of Physical Chemistry B Lipid Order and Dynamics inside the Bilayer. In agreement with lipid packing, there is an increase in order inside the bilayer for the DMPC/cholesterol systems with increasing sterol content (Figures 3A and S3A). Beyond 40%

correlation function C2(t), but this is not the scope of this work as we are interested in the general trends and effects of sterol content and lipid head groups on bilayer properties. It has been well established that high sterol content increases the order of lipids in bilayers,44 as observed by the SCD trends. The sterol molecule also aligns more vertically with respect to the bilayer normal as the order increases, which is exactly what our models show (Figure 4). Additionally, the tilt angle in the

Figure 4. Average sterol chain tilt as a function of cholesterol concentration in both DMPC and DOPC bilayers at 303.15 K. As expected, the sterol aligns more to the bilayer normal as cholesterol concentration increases. The error bars are the standard errors over the three replicates.

DOPC systems is higher than that in the DMPC systems due to the more disordered nature of DOPC compared to DMPC. Table 4 summarizes the cholesterol tilt angles for the other fully saturated systems, DMPA, DMPE, DMPS, DMPG, and DMPI, which is discussed later.

Figure 3. SCD order parameters showing the increase in bilayer order with increasing cholesterol concentration for (A) DMPC and (B) DOPC sn-2 chains at 303.15 K. The standard errors for each measurement are comparable or smaller than the marker size.

Table 4. Average Lipid Tilt Angle in Ld and Solid Gel Bilayers with ± Standard Errors

cholesterol content, there is no significant change in the SCD values. Note that this is the same concentration where the SA levels off, both the overall and component-wise SA (Figure 2). Agreement in these trends was expected; decrease in the SA means the volume that each lipid occupies is also reduced, thus having less freedom to move and increasing the overall order of the tails (higher order parameters). Therefore, with no change in the SA, no significant increase in order is expected in a bilayer from the physical standpoint, which has been previously stated in both experimental and simulation studies.52−54 For DOPC on the other hand, there is no such saturation point in the SCD values, which increase steadily as cholesterol content increases (Figures 3B and S3B). As expected, DOPC is less ordered at all concentrations due to its double bond between C9 and C10 in each tail, which have the lowest order parameter values. The saturation effect on DMPC models is also observed in the lipid wobble (relaxation) times for these systems (Table 3). There is steady increase in average wobble time as cholesterol content increases in DMPC and DOPC membranes. However, there is no statistical difference (within 95% CI) once 40% cholesterol content is reached in the DMPC models, which is not observed for the DOPC models. These DOPC models have constantly increasing average wobble times for increasing sterol content, but there is no statistical difference (within 95% CI) between wobble times, except for the 20% cholesterol simulation. More precise wobble values could be obtained using longer simulation trajectories for the computation of the

system

temp (K)

DMPC DMPA DMPE DMPS DMPG DMPI DMPA2 DMPE2 DMPS2

303.15

333.15

sterol angle (deg) 18.90 15.31 12.78 17.11 20.16 18.55 18.93 17.59 18.08

± ± ± ± ± ± ± ± ±

0.04 0.27 0.63 1.70 1.13 0.01 0.51 0.01 0.02

lipid tail angle (deg) 20.94 13.88 13.61 14.38

± ± ± ±

0.17 0.20 0.51 0.44

22.51 ± 0.29 20.65 ± 0.15 21.13 ± 0.33

A direct comparison of acyl chain order of our systems to experimental data is difficult due to the experimental challenge of individual carbon labeling. As such, ranked order comparison gives a better assessment of the agreement between simulations and experiment. For the DMPC systems with 15% cholesterol, there is fair agreement (within 99% CI for all comparisons) between the calculated values and the experimental ones (Figure 5A), considering the slightly lower sterol concentration in experiment (12.5%). MD simulations of the 20% cholesterol DMPC systems are also in good agreement for the sn-2 chains, yet the sn-1 chains are significantly more disordered in experiment at the high rank, i.e., end of the chain (Figure 5B). This could be due to some phase boundary effects of Lo and Ld domains in the experiment as there is a significant 11765

DOI: 10.1021/acs.jpcb.6b08574 J. Phys. Chem. B 2016, 120, 11761−11772

Article

The Journal of Physical Chemistry B

lipid head groups, albeit less than the changes introduced by altering acyl chain saturation and cholesterol concentration. Although chain order slightly depends on the lipid headgroup, it is clear that the major influence is the chain structure. Lipid Lateral Organization. To further survey the lateral ordering information gleaned from SA per lipid and SCD order parameter data, RDFs were calculated for all systems. When comparing these calculations across replicas and both leaflets of the simulated bilayers, peak heights from the RDFs (Figure S6) are not good metrics due to high levels of variability in molecular clustering and sensitivity to initial conditions (as different initial distributions of lipids on the bilayer plane were used for different replicates). This variability has been outlined in the literature before, suggesting that longer time scales are needed to obtain lateral ordering results that are independent of initial conditions.58 As such, instead of focusing on clustering in individual systems, we focused on trends for the radial location of the first and second solvation shells and the variance in the peaks of these same solvation shells. For DOPC and DMPC, these trends helped to tease out the effects of chain saturation on the solvation shells. The aforementioned high level of variance can be seen in the peak heights for both DMPC and DOPC (Figure 6A,C). However, a pattern can be seen in this variance, with increasing cholesterol concentration leading to a more defined first and second peak height. This effect appears to be less pronounced in DMPC, which may be a result of stronger interactions with the saturated chains, leading to more aberrant clustering behavior. Both DMPC and DOPC simulations display convergence (within 95% CI) to a radial distance of 5.5 Å for the first cholesterol solvation shell, and 11 Å for the second shell (Figure 6B,D). The reduction in the errors of this solvation shell radial positioning as cholesterol concentration is increased is expected, as the low number of cholesterol molecules in the 5% and 15% systems lead to noisier RDFs. The cholesterol-cholesterol RDF statistics for the simulations of the various headgroups with DM chains and cholesterol at 15% are in Figure 7. Along with geometric properties, such as size and shape, the headgroups in the DM systems can be more generally classified as anionic (PA, PG, PI, PS) having a net negative charge, or zwitterionic (PC, PE) having a net neutral charge with distinct positive and negative components within the headgroup. The larger, anionic head groups, PI and PG, show a higher variability in the radial distances of the second solvation shells. This observation appears to be independent of the temperature at which these simulations are run, as the first solvation shell positions show less variability than those of the PC head groups at 15% cholesterol. Figure 8 shows sample snapshots of hydrogen bonding configurations in three specific bilayers; the images show that PS can separate cholesterol using PS−PS hydrogen bonds (Figure 8A) and rarely forms hydrogen bonds with cholesterol (Figure 8B), while PI surrounds the cholesterol molecules and forms a semicircle interconnected by hydrogen bonds (Figure 8D), and PA can form hydrogen bonds between the carbonyl on the acyl chains and cholesterol, albeit rarely (Figure 8C). These exemplary snapshots do not persist through the entire simulation time, although they do persist several nanoseconds and give examples of potentially common interactions between phospholipids and cholesterol. The location of hydrogen bonds varies for lipid−lipid interactions, but lipid−cholesterol interactions occur through hydrogen bonds between the hydroxyl hydrogen of cholesterol and the carbonyl oxygen of

Figure 5. SCD values for DMPC in ranked order, i.e., ordering from the highest value to the lowest SCD regardless of carbon position, showing good agreement between simulation and experiment. Comparisons are grouped by cholesterol concentration (A: 10−15%, B: 20%, and C: 50%). Exp1 at 314.15 K by Martinez et al.;56 Exp2 at 303.15 K by Barry and Gawrisch.57

increase in order at 25% cholesterol (Figure 5B). Temperature is also another factor, with the experimental data of Martinez et al.55,56 (Exp1) being obtained at a temperature of 314 K, which is higher than both simulation and the data of Barry and Gawrisch57 (Exp2) at 303 K. When looking at the raw experimental data of the sn-2 chain in particular, there is a plateau at the top of the chain likely due to the inability to measure a unique signal for these carbons. After this plateau, the experimental data follows the same general trend as our simulations, suggesting that problems with labeling may be present (Figure S4). The comparison with 50% cholesterol content is more favorable (Figure 5C). Both acyl chains for simulation and experiment follow the same pattern in excellent agreement and deviate only slightly at the top ranked chains and near the ends. Comparing the chain order for all carbons as a function of headgroup for DM lipids (Figure S5), the values have only minor differences between head groups and correlate with the SA results. The sn-1 chain shows that the ordering of DMPA and DMPG are similar and differ appreciably from the other 11766

DOI: 10.1021/acs.jpcb.6b08574 J. Phys. Chem. B 2016, 120, 11761−11772

Article

The Journal of Physical Chemistry B

Figure 6. Solvation shell characterizations from RDFs for cholesterol-cholesterol highlight the complex sterol interactions in (A,B) DMPC and (C,D) DOPC simulations in terms of solvation shell peak height (left) and solvation shell radial position (right).

Figure 7. Solvation shell characterizations from RDFs for cholesterol-cholesterol highlight the complex sterol interactions in head groups with saturated chains in the Ld phase (303.15 K for DMPI and DMPG; 333.15 K for DMPA, DMPE, and DMPS) in terms of solvation shell peak height (A) and solvation shell radial position (B).

the lipid. Table 5 summarizes the total number of hydrogen bonds in all DM systems. Due to the relatively low concentration of cholesterol in the bilayer, lipid−cholesterol interactions are relatively infrequent, but their occurrences coupled with lipid−lipid interactions are possibly what stabilizes the secondary solvation shells of the RDFs. Lipid Membrane Profile Densities. Head-to-head thickness (DHH), overall bilayer thickness (DB), and hydrophobic thickness (2DC) from simulation data can all be determined using the EDPs and volume probability plots for the different lipid sections (tails, head groups, phosphate, and methyl regions, Figure S7). Fourier transforms of the simulation data used to calculate the EDPs can be directly compared to X-ray form factors from experiment using SIMtoEXP software.42 Using this data, DMPC shows excellent agreement at 5% and 20% cholesterol (p-value >0.95, χ2 goodness-of-fit test) and almost no agreement at 40% cholesterol (p-value