Improvement of Parameters of the AMBER Potential Force Field for

Jun 24, 2015 - In this study, we improved parameters of the AMBER potential force field for phospholipids in order to describe the thermal phase trans...
0 downloads 4 Views 9MB Size
Article pubs.acs.org/JPCB

Improvement of Parameters of the AMBER Potential Force Field for Phospholipids for Description of Thermal Phase Transitions Koji Ogata* and Shinichiro Nakamura Nakamura Laboratory, RIKEN Innovation Center, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan S Supporting Information *

ABSTRACT: In this study, we improved parameters of the AMBER potential force field for phospholipids in order to describe the thermal phase transition using molecular dynamic (MD) simulations. To estimate the errors of the main phase transition temperature (Tm), first, MD simulations using the GAFFlipid and Gaff parameters were performed for six phospholipid bilayers, 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC), 1,2-dipalmitoyl-sn-glycero-3-phosphoethanolamine (DPPE), 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC), 1,2-distearoyl-sn-glycero-3-phosphocholine (DSPC), 1-palmitoyl,2-oleoyl-sn-glycero-3-phosphocholine (POPC), and 1-palmitoyl,2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE), with increasing temperature. The Tm values were characterized according to the structural parameter, area per lipid, and gauche ratio in alkyl chains. The Tm values of the six lipids showed ∼50 K differences from the experimentally measured values. To reduce these errors, the well-depth values in the Lennard−Jones potential of the alkyl chains were modified to fit the Tm values of the simulation to the experimental values in a single DPPC bilayer. After the fitting procedure, the Tm values of the six lipids improved, and the errors of Tm improved from ∼50 to ∼15 K. We show that the simulation applying the improved parameters provides more accurate results than the original parameters. These modified parameters were also found to be useful for performing MD simulation of transmembrane proteins with membrane models.

1. INTRODUCTION The plasma membrane, which is composed of a lipid bilayer mainly consisting of phospholipids, sphingolipids, and cholesterol, plays an important role in governing cellular processes, such as adhering cells, conducting ions, and transforming signals by transmembrane proteins.1 Understanding the physicochemical properties of the plasma membrane is of critical importance in the chemical, biological, and medicinal science fields. The thermal behaviors of lipids are investigated using various methods, e.g., differential scanning calorimetric and X-ray diffraction studies.2−4 Many scientists have studied the phase transition of a single lipid bilayer using these methods. For instance, a single DPPC bilayer has four phases, crystalline (Lc), gel (Lβ), ripple (Pβ), and liquid crystalline (Lα), and the temperature of this phase transition is well known to be ∼20 °C (∼293.15 K) for Lc → Lβ (subphase transition), ∼35 °C (∼308.15 K) for Lβ → Pβ (prephase transition), and ∼41 °C (∼314.15 K) for Pβ → Lα (main phase transition).2,3,5 For each phase, simulations applying experimental data (NMR, FTIR, neutron scattering, etc.) have provided important information in various fields.6−8 However, the parameters still need improvements to reproduce the thermal phase transition consistently, due to concerns of the phase transition temperature changing in the simulation. © XXXX American Chemical Society

Simulations on lipid bilayers consisting of single lipids, especially DPPC and DPPE, have been studied using coursegrain models,9,10 united atom,11 and all atoms.12 These studies have been reported to review the thermodynamic behavior of the lipid bilayer. Also, in our previous study, MD simulation using the AMBER potential force field successfully linked the experimentally obtained information on the structures of single DPPC and DPPE bilayers in each phase to all of the physicochemical properties.13 We reported the mechanism of the gel and ripple phases formations in the DPPC bilayer. Our results indicate that the force generated by the conformational changes in the alkyl chains in the Lβ phase is essential for forming ripple conformation in the Pβ phase. However, in these simulations, the error of the Tm was ∼50 K, compared with the experimentally measured values. This error raises significant problems when performing simulations of the lipid bilayer with transmembrane proteins. In this study, we improved parameters of the AMBER force field for lipids to reproduce the Tm value in MD simulation. Previous studies on the improvement of the parameters of the force field focused on the lipid conformation at the Lβ or Lα Received: February 18, 2015 Revised: June 23, 2015

A

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B phase.14−18 Lipid conformations using these parameters showed good agreement with the experimental conformations; however, these parameters did not reproduce the Tm value. It is generally known that MD simulations overestimate the Tm values. Therefore, it was not clear whether the phase identified from the lipid conformation corresponds to the given temperature. To avoid this problem, the force field parameter sets are fitted to reproduce the Tm value as accurately as possible. In particular, in the case of the Gaff parameter set, the parameters need to be improved. In this study, the well-depth values in Lennard−Jones potential in alkyl chains were modified to express the Tm value as closely as possible to the experimental ones. We validated the improved parameters using MD simulation while increasing the temperature. The results demonstrate that the error improved to within 15 K, and the thermal behavior in each phase was similar to that in our previous study. We applied the modified parameters to several single lipid bilayers, DOPC, POPC, DSPC, DPPE, and DOPE. The results for these six lipids reproduced Tm values consistent with the experimental data.

snapshots were obtained in total. The 2000 snapshots were classified by a clustering method with a 1.0 Å threshold value. In each cluster, the conformation with the lowest energy was regarded as the representative conformation. Finally, 1130, 1235, 1026, 1107, 1436, and 1026 conformations were obtained for DPPC, DSPC, POPC, DOPC, DPPE and POPE, respectively. ESP charges were obtained by single-point calculation (by Gaussian0923) at the B3LYP/6-31G(d) level with the polarizable continuum model (PCM). The partial charge for the nth atom is given as cn =

GAFFlipid

a

z

area per lipid (Å2)

no. of waters/lipid

DPPC DSPC POPC DOPC DPPE POPE

64 64 64 68 64 64

64 64 66 68 64 66

106 107 125 96 103 97

64 64 66 72.25 64 66

35.32 37.13 37.55 37.22 37.51 37.39

Gaff

atom type

r0 (Å)a

ε (kcal/mol)b

r0 (Å)a

ε (kcal/mol)b

alkyl c3 alkyl hc glycerol c3 glycerol hc

2.010 1.340 1.908 1.487

0.055 0.024 0.1094 0.0157

1.908 1.487 1.908 1.487

0.1094 0.0157 0.1094 0.0157

Van der Waals radii. bWell depth.

parameter fitting procedure for both vdW parameters. Note that the differences between GAFFlipid and Gaff are only in vdW parameters in the alkyl chains in this study. 2.2. MD Simulation. We performed MD simulation for six lipids in three stages. First, the Tm values by the original parameters were measured. Then, short period simulations for parameter fitting were performed. Finally, simulations were performed with the fitted parameters. The Tm values were evaluated by the fitted and original parameters. The lipid bilayer models were optimized using the conjugate gradient method in 10 000 steps for all atoms after 10 000 steps for hydrogen atoms using AMBER12.2,24,25 After optimization, 100 ps MD simulation under constant volume (NVT) was performed with constraint for the lipids. In this study, all MD simulations were carried out with the 10 Å cutoff and particle mesh Ewald (PME) method26,27 in the periodic boundary condition. To generate an initial structure, 10 ns MD simulation under constant pressure (1 atm) and temperature (NPT) using a Langevin thermostat28 with a collision frequency of γ = 1 and no surface tension was performed. In this procedure, anisotropic Berendsen control29 was used with a time constant of 1 ps. The MD simulations for the six lipid bilayers were started at the temperature, as shown in Table 3. The simulations were performed for 70 (for parameter fitting procedure), 320 (for original parameter), or 260 ns (for fitted parameter) while increasing the temperature 5 K every 10 ns. The snapshots

system size (Å) y

(1)

i=1

Table 2. Lennard−Jones Parameters in GAFFlipid and Gaff

Table 1. Initial Structure of Lipid Bilayers x

N

∑ cin exp{−Ei /kT }

where ci and Ei are the partial charge and energy at the ith conformation, respectively, k is Boltzmann’s constant, and T is temperature. Z is summed over exp{−Ei/kT}. These parametrizations of the lipids were applied to DPPC, DSPC, POPC, DOPC, DPPE, and POPE (Figure S1, Figure S2, Table S1, and Table S2, Supporting Information). The parameters for vdW interaction in the lipid headgroup were obtained from Gaff. In the alkyl chain, the carbon and hydrogen atoms in GAFFlipid and Gaff were slightly different, as shown in Table 2. Therefore, we tested and performed a

2. MATERIALS AND METHODS 2.1. Modeling of Lipid Bilayer and Parameter of Lipids. Lipid bilayer models were built using the same protocol as our previous study.13 The lipid models were constructed with extended alkyl chains in which the torsion angles of saturated and unsaturated bonds were set at 180° and 0°, respectively, so that the chains of both form a parallel shape. These lipid models were aligned along the z axis and the bilayer on the xy plain. The bilayer models contain 128 lipids (8 × 8 × 2). The area per lipid of the initial model was set similar to the experimental data measured in the Lα phase, as shown in Table 1. The lipid bilayer models were hydrated into the box of TIP3P water models.19

name of lipids

1 Z

The atom types in GAFF20 were assigned to the lipids using an antechamber software package.21 The parameters for the torsion angle energy in the alkyl chains were obtained from GAFFlipid reported by Dickson et al.22 This parameter set was well improved for the alkyl chain in the lipids. The other parameter values for bonded energies, i.e., bonds, angles, and torsion angles in the head groups and improper torsion angles, were obtained from Gaff. In GAFFlipid, the DSPC and DPPE parameters are not available to the public. To maintain the consistency of all parameter sets, we calculated the atomic charges in the same way as our previous study. First, the initial partial charges were obtained by using an antechamber at the am1bcc level. The lipids were relaxed by optimization in the vacuum state. After optimization, unconstrained MD simulation at a temperature of 300 K was performed for 2 ns, and snapshots were stored at time intervals of 1 ps. Two thousand B

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

case of Gaff, we tested all combinations of (αc3, αhc), in which the αc3 and αhc values were varied in 4/8, 5/8, 6/8, 7/8, and 8/ 8. To find the optimized α value, the first screening with short period MD simulation was performed for various αc3 and αhc values. The temperature in the simulation increased from 303 to 328 K. Using optimal well-depth values, the long period MD simulation was performed for the six lipid bilayers, and the thermal behavior of these lipids was analyzed. 2.4. MD Simulation Using LIPID14, CHARMM, and GROMACS. To compare the currently modified parameters with other force fields, MD simulations with lipid14,30 CHARMM,31 and GROMACS32 were performed for the DPPC and DPPE single lipid bilayer. In the case of simulation using lipid14, the same models and the same procedure, as mentioned above, were used. In the simulation of CHARMM, the chamber software,33 contained in the AMBER software package, was used with the CHARMM36 force field for lipids.17 The simulation was performed under the same conditions as our work. In the simulation using GROMACS, the GROMACS native engine32 was used with the GROMOS 53a6 force field15 and with the modified vdW radius in the lipid headgroup.34 In the solvation, the SPC water model35 was used instead of TIP3P. The Parrinello−Rahman barostat36 was applied instead of Berendsen control. 2.5. Order Parameters. To validate the parameter performance, the profile of the deuterium-order parameter of the alkyl chain, SCD, which is a measure of the orientation of the lipids with respect to the bilayer normal, was evaluated from the MD trajectories. Here, SCD is calculated as

Table 3. Temperature Ranges in MD Simulations temperature (K) original parameters

fitted parameters

name of lipids

start

end

start

end

DPPC DSPC POPC DOPC DPPE POPE

263 263 223 203 263 223

418 418 378 358 418 378

243 243 203 183 243 223

368 368 328 308 368 348

were stored every 10 ps. The temperature dependency of conformational behavior was monitored by snapshots. 2.3. Parameter Fitting Procedure. We performed the parameter fitting procedure to improve the Tm of a single DPPC bilayer in the simulation. The well-depth values in Lennard−Jones potential in the alkyl chains were modified so that the Tm of the simulation matches the experimental ones (∼314.15 K). Here, the nonbonded interaction consisting of the vdW and coulomb interaction is given by U=

qiqj εrij

+

⎧⎛ 0 ⎞12 ⎛ r 0 ⎞6 ⎫ ⎪ rij ⎪ ij ⎜ ⎟ eiej ⎨⎜ ⎟ − 2⎜⎜ ⎟⎟ ⎬ ⎪⎝ rij ⎠ ⎝ rij ⎠ ⎪ ⎩ ⎭

(2)

0

where rij and rij are distance and standard vdW radius, respectively, between atom i and atom j. ε is the dielectric constant, and qi and ei are partial charges and well depth at atom i, respectively. The well depth of atoms i and j represents the minimum energy of the vdW interaction between atoms i and j. Here, the ei values of the sp3 carbon and the hydrogen atom in the alkyl chains were changed (Table 2). In an attempt to find appropriate scaling factors multiplying those well-depth values, the following track was taken. Here the modified well depths are presented as ec3′ = αc3·ec3 and ehc′ = αhc·ehc. Then, the scaling factor pair components c3 and hc are denoted as (αc3, αhc). In the case of GAFFlipid, the ec3 value is already small. Therefore, only ehc was scaled, and attempts were made to find the appropriate αhc (αhc < 1) . In this study, 8 cases, αhc ={0/8, 1/8, 2/8 3/8, 4/8, 5/8, 6/8, 7/8}, were tested. In the

SCD =

1 ⟨3 cos2 θCD − 1⟩ 2

(2)

where θCD is the angle between a vector of the C−H bond and the bilayer normal and ⟨⟩ denotes an ensemble average.37,38 The SCD values were calculated and compared with the experimentally measured data.

3. RESULTS 3.1. Tm Values of Six Lipids Using Original Parameters. 3.1.1. Tm Values of Six Single Lipids Using Original GAFFlipid.

Figure 1. Temperature vs (a) area per lipid (AL), (b) gauche ratio in fatty acid (RF), (c) differential of AL (|ΔAL|), and (d) differential of RF (|ΔRF|) of DPPC, DOPC, POPC, DSPC, DPPE, and POPE using GAFFlipid. C

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

however, were in agreement with the experimental values. For instance, the difference between DPPC and DSPC in the simulation was ∼15 K, which was in good agreement with the experimental relative value (∼14 K). The value between DSPC and DPPE in the simulation was ∼15 K. This value was slightly larger than the experimental value (∼8 K) but within the acceptable difference compared to the Tm value differences. As a whole, the absolute values of the Tm showed a difference of ∼50 K in the six lipids, but the relative values of the Tm values more or less agreed with the experimental data. 3.1.2. Tm Values of Six Single Lipids Using Gaff. MD simulation was carried out to measure the error of the original parameter of Gaff. The overall profiles of the AL and RF for the six lipids obtained using Gaff were completely different from those using GAFFlipid (Figure 2). The Tm values of Gaff were much larger than those of GAFFlipid. The main phase transition for the PC headgroup was detected in the same temperature range with GAFFlipid but not for the PE headgroup. The Tm value of the PC headgroup showed a large error, ∼90 K different from the experimental value (Table 4), while that of the PE headgroup was over 80 K, although the parameter of the vdW interaction in the headgroup and charges in the coulomb interaction were the same values for GAFFlipid and Gaff. Therefore, the vdW interaction in the tail group of the lipids may cause the difference in the Tm values. Comparison of the vdW parameters between GAFFlipid and Gaff showed that the differences in vdW radius and well-depth values were slight, but the thermal behavior was significantly different in the simulations. These small differences therefore were very effective and sensitive to the results of simulations. The torsion angle parameters used in GAFFlipid may have contributed to the improvement of the results of the simulation for lipid bilayers. In other words, the original torsion angles in Gaff may be one of the causes of the errors. 3.2. Results of Parameter Fitting. 3.2.1. Parameter Fitting for GAFFlipid. The parameter fitting procedure for GAFFlipid was successfully performed using the short period MD simulation with increasing temperature. In some cases of the simulations, the profiles of AL and RF showed phase transition (Figure 3). In Figure 3, the entire profile of the simulation using the original parameter values, (αc3, αhc,) = (1,

In order to measure the error of the original parameters of GAFFlipid, we performed MD simulations. In the MD simulations using GAFFlipid, the phase transition of the six lipids was observed. The phase transition is traced by area per lipid (AL), gauche ratio in fatty acid (RF), |ΔAL| (= |AL(T) − AL(T − 5)|; T, temperature), and |ΔRF| (= |RF(T) − RF(T − 5)|), as shown in Figure 1. In particular, DPPC, DPPE, and DSPC, which have two palmitic acids (16:0) or two stearic acids (18:0), showed a clear phase transition at 358, 388, and 373 K, respectively. Also, POPC and POPE, which consists of one palmitic acid and one oleic acid (18:1), showed phase transition at 313 and 343 K, respectively. On the other hand, DOPC, having two oleic acids, showed a gradual increase of AL and RF, which is not the profile of clear phase transition. When the profiles of POPC and POPE are compared, the phase transition in POPE is more distinct than POPC. The Tm value of POPE was ∼20 K higher than that of POPC. It is expected that, in the simulation, the higher the Tm of the lipids, the clearer the phase transition behavior of the lipid. The Tm values of the six lipids in the simulation shifted at a higher temperature range than the experimentally measured values (Table 4). The errors of the Tm values for DPPC, DSPC, Table 4. Summary of the Main Phase Transition Temperature experimenta

a

GAFFlipid

Gaff

name of lipids

Tm (°C)

Tm (K)

Tm (°C)

Tm (K)

Tm (°C)

Tm (K)

DPPC DOPC DSPC POPC DPPE POPE

41 −20 54 −2 63 25

314.15 253.15 328.15 271.15 336.15 298.15

84.85 29.85 99.85 39.85 114.85 69.85

358 303 373 313 388 343

129.85 44.85 144.85 79.85 144.85< 104.85
418 >378

Tm values are approximated values.

and DPPE were ∼43, ∼45, and ∼54 K, respectively. Also, those of POPC and POPE were ∼40 and ∼45 K, respectively. These errors were significantly large. The relative values in the Tm,

Figure 2. Temperature vs (a) area per lipid (AL), (b) gauche ratio in fatty acid (RF), (c) differential of AL (|ΔAL|), and (d) differential of RF (|ΔRF|) of DPPC, DOPC, POPC, DSPC, DPPE, and POPE using GAFF. D

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 3. Temperature vs (a) area per lipid (AL), (b) gauche ratio in fatty acid (RF), (c) differential of AL (|ΔAL|), and (d) differential of RF (|ΔRF|) of DPPC using modified well-depth parameters in GAFFlipid. Simulations were performed in the temperature range from 303 to 323 K.

8/8), is also shown as reference. The Tm value of (1, 8/8) was at 358 K, which is ∼43 K larger than that of the experimentally measured one, and the AL and RF values at Tm were 56.8 Å2 and 0.24, respectively. As shown in Figure 3, the αhc values strongly influence the Tm value in the simulations. In the range of αhc ≤ 3/8, AL and RF were close to the values in the Lα phase obtained by original parameters. The surface of the lipid bilayer was expanded, and the alkyl chains were completely melted. In this range, the main phase transition might appear below 303 K. The small value of αhc therefore results in the low Tm value. In this case, the alkyl chains easily form the gauche formation. Consequently, the DPPC bilayer does not show phase transition near 313 K. On the other hand, in the range of αhc ≥ 5/8, AL and RF showed relatively small values, which were ∼45 Å2 and ∼1, respectively, at 303 K. The AL and RF of αhc ≥ 5/8 were very similar to those of αhc= 8/8 in the low-temperature range. According to our previous study, the DPPC bilayer corresponding with the AL and RF values was in the Lc phase.13 Therefore, in the range of αhc ≥ 5/8, ehc′ was too large to overcome the barrier to take the gauche formation. Consequently, the lipid bilayers still remained in the Lc phase. The case of αhc = 4/8 (ehc′ = 0.012) was close to the experimentally measured result. The values of AL and RF at 318 K were 57 Å2 and 0.24, respectively. Those values were similar to the original parameter values (56.8 Å2 and 0.24) at 358 K, which correspond to the original parameters of the Tm in the simulation. Therefore, the long period MD simulations were performed using αhc = 4/8 for the six lipids. The modified GAFFlipid is designated as mod_GAFFlipid in the following sections. 3.2.2. Parameter Fitting for Gaff. Figure 4 shows a summary of the Tm values for the different scaling factors of (αc3, αhc) pairs using Gaff. Figure 5 shows the profiles of AL and RF with (αc3, αhc) pairs, which indicate that the Tm value is in or near the temperature range of the fitting procedure. The Tm value was very sensitive to the pair of (αc3, αhc) values (Figure 4). A small change in the (αc3, αhc) pair showed dramatically different Tm values. For instance, in comparing the results of (5/8, 5/8) and (6/8, 5/8), both Tm values did not appear in the temperature range. In the simulation of (5/8, 5/8), the bilayer surface had already started to expand and melt the alkyl chains

Figure 4. Summary of Tm values on scaling factors of (αc3, αhc) in Gaff. Blue and red circles indicate that the Tm value is below 303 K and over 328 K, respectively. Numbers in the ellipse are the Tm values when the scaling factor is used.

at 303 K, whereas in the simulation of (6/8, 5/8), the bilayer surface was starting to expand and melt the alkyl chains at 328 K (Figure 5). On the basis of these results, we selected two pairs, (6/8, 4/ 8) (ec3′ = 0.08205 and ehc′ = 0.00785) and (4/8, 8/8) (ec3′ = 0.0547 and ehc′ = 0.0157). Here those pairs were called mod_Gaff and mod_Gaff2, respectively. The well-depth values in these pairs are opposite to each other, i.e., mod_Gaff has a large ec3′ and a small ehc′, whereas mod_Gaff2 has a small ec3′ and a large ehc′. For the two pairs, the long period MD simulations with increasing temperature were performed. The effects of these differences on the thermal behaviors in the lipids were observed in our simulation, as mentioned below. 3.3. Simulation Using Modified Parameters. 3.3.1. Tm Values of Six Lipids Using Modified Parameters. The entire profiles of AL and RF of the three modified parameters obviously shifted to the lower temperature range (Figures 6, 7, and 8) when compared with Figure 1. The peak temperature of DPPC in the graphs on the three parameters was found to be between 303 and 323 K, which is the same range as the fitting procedure (vide supra). Therefore, the short period MD simulation is useful for find the phase transition temperature. E

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 5. Temperature vs (a) area per lipid (AL), (b) gauche ratio in fatty acid (RF), (c) differential of AL (|ΔAL|), and (d) differential of RF (|ΔRF|) of DPPC using modified well-depth parameters in Gaff. Simulations were performed in the temperature range from 303 to 323 K.

Figure 6. Temperature vs (a) area per lipid (AL), (b) gauche ratio in fatty acid (RF), (c) differential of AL (|ΔAL|), and (d) differential of RF (|ΔRF|) of DPPC, DOPC, POPC, DSPC, DPPE, and POPE using mod_GAFFlipid.

Figure 7. Temperature vs (a) area per lipid (AL), (b) gauche ratio in fatty acid (RF), (c) differential of AL (|ΔAL|), and (c) differential of RF (|ΔRF|) of DPPC, DOPC, POPC, DSPC, DPPE, and POPE using mod_Gaff.

F

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 8. Temperature vs (a) area per lipid (AL), (b) gauche ratio in fatty acid (RF), (c) differential of AL (|ΔAL|), and (d) differential of RF (|ΔRF|) of DPPC, DOPC, POPC, DSPC, DPPE, and POPE using mod_Gaff2.

Table 5. Summary of the Main Phase Transition Temperature Using Modified Parameters experimenta

a

mod_GAFFlipid

mod_Gaff

mod_Gaff2

name of lipids

Tm (°C)

Tm (K)

Tm (°C)

Tm (K)

Tm (°C)

Tm (K)

Tm (°C)

Tm (K)

DPPC DOPC DSPC POPC DPPE POPE

41 −20 55 −2 63 25

314.15 253.15 328.15 271.15 336.15 298.15

49.85 −30.15 54.85 9.85 64.85 29.85

323 243 328 283 338 303

44.85 −20.15 49.85 14.85 64.85 9.85

318 253 323 288 338 283

39.85 −30.15 49.85 19.85 54.85 19.85

313 243 323 293 328 293

Tm values are approximated values from experimental data.

Table 6. AL in Lα Phase for Six Lipids experiment

a

MD simulation

name of lipids

T (°C)

T (K)

AL (Å2)

T (K)

mod_GAFFlipid, AL (Å2)

mod_Gaff, AL (Å2)

mod_Gaff2, AL (Å2)

DPPC DOPC DSPC POPC DPPE POPE

50 30 65 30 69 30

323.15 303.15 338.15 303.15 342.15 303.15

63.3a 67.4,b 72.5c 66a 64.3,d 68.3e 60.5a 59f

323 303 338 303 343 303

56.2 63.4 60.9 58.8 56.4 52.1

61.2 63.3 62.6 60.4 58.5 55.8

63.1 65.7 65.5 61.9 62.1 58.9

Petrache et al.39 bKučerka et al.6 cNagle and Tristram-Nagle.4 dKučerka et al.40 eKučerka et al.41 fRappolt et al.42

values; their errors are within ∼15 K. This value was relatively smaller than that of the original GAFFlipid (∼40 K) and Gaff (over 80 K). The largest error was ∼20 K, which was observed in POPC using mod_Gaff2, still smaller than that by the original parameter (∼80 K). In this case, the phase transition was unclear because two different temperatures were obtained; | AL| showed the peak at 293 K, whereas |RF| was at 278 K. Therefore, from the AL and RF profiles, 293 and at 278 K are candidates of the Tm of POPC. For the three parameters, the profiles of AL and RF for DOPC and POPC were found to gradually increase. It is hard to identify the main phase transition. Small peaks were also found from |ΔAL| and |ΔRF|. The experimentally measured Tm values of these lipids are below 0 °C (278.15 K). In this temperature range, the water molecules have small motion. Therefore, in the low-temperature range, the energy transferred from the thermal vibrations seems not to accumulate enough to

Note that the 100 ns MD simulations, which were additionally performed for 90 ns to the 10 ns simulations of DPPC at each temperature, showed the same characteristics as those of 10 ns MD simulations (Figure S3, Supporting Information). The Tm values in the 100 ns MD simulations are with ∼10 K difference from those of the 10 ns MD simulations. The profiles at the low- and high-temperature ranges, however, did not differ significantly. The order parameters in the 100 ns simulations were also very similar to those of the 10 ns simulation (Figure S4, Supporting Information). These results show that the continuously performed longer period simulation at each temperature produces the same results as this work. Table 5 shows the Tm values identified from the AL and RF profiles. The Tm values in all simulations were distinctly improved. They are much closer to the experimentally measured values than those of the original parameters (Table 4). Almost all simulations show the improvement of the Tm G

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 9. Lateral radial distribution function, g(r), between lipids in simulation using mod_Gaff2. Profiles were drawn in gradational color from black to gray corresponding to the temperature from low to high.

Figure 10. Order parameter of (a) sn-1 and (b) sn-2 in DPPC, (c) sn-1 and (d) sn-2 in POPC, and (e) sn-1 and (f) sn-2 in DSPC. The values of Exp_1, Exp_2, and Exp_3 were taken from reports by Leftin and Brown,43 Douliez et al.,44 and Leftin et al.,45 respectively.

mod_Gaff2 were 56.2, 61.2, and 63.1 Å2, respectively, at 323 K. The experimental values of AL are 63.3 Å2 at 323.15 K. The values of mod_Gaff2 were closer to the experimental one than those of mod_GAFFlipid and mod_Gaff. This reasonable feature by mod_Gaff2 can be found for all lipids as shown in Table 6. Comparing the well-depth values for mod_Gaff and mod_Gaff2, the ec3′ of mod_Gaff is larger than that of mod_Gaff2, whereas ehc′ of mod_Gaff is smaller than that of mod_Gaff2. For both mod_GAFFlipid and mod_Gaff, the total sizes of the vdW radius in carbon and hydrogen atoms were

expand the bilayer surface near the Tm. Instead, the energy may transfer to the bilayers. As a consequence, the lipid surface gradually expanded and the phase transition temperature cannot be clearly identified. The AL and RF profiles were very similar for the three modified parameter sets. By contrast, remarkably different points were found in the Lα phase after Tm (Table 6). In the Lα phase, the increase in AL by mod_GAFFlipid and mod_Gaff is somewhat smaller than those by mod_Gaff2. For instance, the AL values of DPPC with mod_GAFFlipid, mod_Gaff, and H

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 11. Electron density profile of (a) DPPC (323 K), (b) DSPC (338 K), (c) POPC (303 K), (d) DOPC (303 K), (e) DPPE(338 K), and (f) POPE (313 K).

their values are close to the experimentally measured ones without significant errors in bilayer conformations such as the packing of lipids. 3.3.3. Analysis of Order Parameters. The SCD profiles of the saturated alkyl chains demonstrate similar characteristics for the three improved parameters, as shown in Figure 10. Still, the SCD profiles of POPC at sn-2 (unsaturated alkyl chain) can show obvious difference for the three parameters. The carbon atoms from C2 to C8 in mod_GAFFlipid show the largest values for the other parameters, whereas those from C12 to C18 are the smallest. Since the SCD profiles at sn-1 (Figure 10c) are very similar for the three parameters, the alkyl chains including C C bonds strongly contribute to make these differences. According to the RF profiles, POPC gradually melts and did not easily show the main phase transition (Figure 6b). In particular, POPC with mod_GAFFlipid did not show the main phase transition clearly in the RF profile. The large SCD value means that the alkyl chain is perpendicular for the xy plane. These alkyl chains were often found from the extended lipids in the Lc and Lβ phases. The simulation of POPC with mod_GAFFlipid therefore has a larger number of extended lipids than those with the other parameters. The profiles of SCD for the three parameters show similar characteristics with the experimentally measured ones. Meanwhile, each value of these profiles was slightly larger than the experimental ones. In particular, the carbon atoms near the carbonyl moiety, from C2 to C9, tend to be overestimated in the saturated alkyl chains. These results imply that the slightly larger number of alkyl chains is still extended at the parts near the carbonyl moiety. These results were independent of the simulation times, since similar profiles were observed from the 100 ns MD simulations. Then the SCD profiles obtained at 15 K higher temperatures than those in Figure 10 were compared with the experimental values. These profiles showed agreement with the experimentally measured data (Figure S7, Supporting

identical. The vdW radius of the carbon for GAFFlipid is larger than that of mod_Gaff, whereas the vdW radius of the hydrogen in GAFFlipid was smaller than that of mod_Gaff. The larger well depth and vdW radius indicate the larger contribution of the thermal behaviors. These results confirm that the parameter for the hydrogen atom having large values of vdW radius and well depth provides a conformation that is more natural in lipid bilayers in the Lα phase. 3.3.2. Lateral Radial Distribution Function. The lateral radial distribution function, g(r), showed clear peaks in the lowtemperature range and broad peaks in the high-temperature range (Figure 9). In the low-temperature range, the lipid bilayer was in the Lc and/or Lβ phases. These phases contained a large number of lipids with extended alkyl chains, which helped maintain the positions of each lipid during the simulation and maintain relative distances between the lipids. On the other hand, in the high-temperature range, most of the alkyl chains of the lipids melted and the lipids defused on the bilayer. Therefore, at distances more than 10 Å, the g(r) profile was constant. The lipids were chaotically distributed in the bilayer. Overall, the g(r) profiles with modified parameters were close to those with the original parameters, GAFFlipid and Gaff (Figure 9, Figure S5 and Figure S6, Supporting Information). The positons of the peaks were very similar among five parameters. The height of the peak of Gaff, however, slightly differed from the other parameters. Simulation using Gaff was conducted in the Lc and Lβ phases at most temperatures, even in the high-temperature range. In this case, the lipids moved slightly in the bilayer and kept their relative distances. Therefore, the profiles of Gaff showed clearer peaks than those of other parameters. On the other hand, the packing of the lipids in the bilayer was similar for all parameters. These results reveal that the packing conformations of the lipids in the simulation using the modified parameters were similar to those using original parameters. Therefore, the modified parameters with rescaled well-depth values can improve the Tm value, as I

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 12. Electron density profile of DPPC using (a) mod_GAFFlipid, (b) mod_Gaff, and (c) mod_Gaff2. Electron density profiles of water and fragments in DPPC, which are the same as those reported by Kučerka et al.,40 are also shown. Methyl groups in choline are described as CholCH3, and the phosphate and part of the choline (CH2CH2N) is P. Carbonyl and glycerol groups is CG, and hydrocarbon in alkyl chains are CH2 and CH3, respectively.

Figure 13. Density profile for water molecules of (a) DOPC, (b) DPPC, (c) DPPE, (d) DSPC, (e) POPC, and (f) POPE using mod_Gaff2. Profiles were drawn in gradational color from black to gray corresponding to the temperature from low to high.

our simulations using the improved parameters are able to reproduce the original state of the lipid bilayers. The fragment-based electron density profiles of DPPC were almost similar for the three modified parameters, as shown in Figure 12. The electron density profiles of DSPC, POPC, and DOPC also showed similar characteristics for the three parameters (Figure S8, Supporting Information). Comparing these profiles with the experimental ones reported by Kučerka et al.,6,40,41 who obtained the electron density profiles from the models by simulations combining X-ray and neutron scattering data, the profiles using the modified parameters are similar to the experimental ones. The peaks of CG, P, and CholCH3 in DPPC of the simulations were around −17, −21, −21 Å, respectively, whereas those of the experimentally measured data were approximately around −14, −20, and −20 Å, respectively. These results reveal that the positions of the fragments in the simulation were slightly distant from the center of the lipid bilayer than the experimental ones. On the profiles of the other lipids, the fragments also tend to be located further from the center of the bilayer (Figure S8, Supporting Information).

Information). Therefore, the alkyl chains have not melted enough in the direction vertical to the lipid bilayer. 3.3.4. Analysis of the Electron Density Profile. The total electron density profile was calculated from the MD trajectories for the six lipid bilayers, as shown in Figure 11. The profiles with the three improved parameters showed a slightly different feature. The profile with mod_GAFFlipid showed sharp features near the peaks, and on the whole, the electron density is higher than those of the other parameters, whereas the overall electron density in mod_Gaff2 is slightly lower than those of the other parameters. The overall shapes of the electron density profiles with the simulations were similar to the experimental findings.41,46,47 The distance between peaks, DHH, of DPPC, DOPC, POPC, and POPE were ∼40, ∼39, ∼39, and ∼40 Å, respectively. These values are slightly larger than those of experimentally measured values of which DPPC, DOPC, POPC, and POPE were 38.3,4 36.9,4 37,41 and 39.5 Å,42 respectively. We believe that these are within the acceptable differences. For this reason, J

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 14. Temperature vs (a) area per lipid (AL), (b) gauche ratio in fatty acid (RF), (c) differential of AL (|ΔAL|), and (d) differential of RF (|ΔRF|) of DPPC using mod_Gaff2, lipid14, CHARMM, and GROMACS.

Figure 15. Temperature vs (a) area per lipid (AL), (b) gauche ratio in fatty acid (RF), (c) differential of AL (|ΔAL|), and (d) differential of RF (|ΔRF|) of DPPE using mod_Gaff2, lipid14, CHARMM, and GROMACS.

at Tm. At a distance of around ±20 Å, small peaks were observed in most of the profiles. By analyzing the fragmentbased electron density profile (Figure 12 and Figure S8, Supporting Information), the phosphoric acid group (P) was located at a similar position. Therefore, these peaks were observed by the water molecules interacting with the phosphoric acids. These interacting water molecules help realize stable bilayer formations during the simulation.

These differences were due to the various distances between peaks. 3.3.5. Analysis of Density Profile of the Water Molecule. To analyze the water−lipid interaction, the density profiles of water molecules were analyzed. Overall, the density profile of the water molecule in the GAFFlipid was similar to the three modified parameters (Figure 13, Figure S9 and Figure S10, Supporting Information). However, the profile of Gaff was slightly different from the other parameters at (αc3, αhc) = (1, 8/8). The simulation using Gaff showed the Tm value at high temperatures. Some of them remained in the Lc and Lβ phases and did not exhibit phase transition in the simulation. Therefore, the lipid bilayer was extended during the simulation at most temperatures and showed a slightly different profile. In Figure 13, the profiles of water molecules of the simulation are obviously separated into two groups depending on the temperature. These separations can be found in the profiles with the different parameters (Figures S9 and S10, Supporting Information). These groups, appearing at lower and higher temperature separately, could be assigned to the series of the Lc and/or Lβ phases and series of the Lα phase, respectively,

4. DISCUSSION To compare the lipid conformations by other potential force fields, MD simulations using lipid14,22 CHARMM,31 and GROMACS32 were performed for single DPPC or DPPE bilayers. The main phase transitions of DPPC and DPPE bilayers were clearly identified from the AL and RF graphs in the simulations using four potential force fields (Figures 14 and 15). The simulation of DPPC using mod_Gaff2 indicated the Tm values near the experimental Tm value (314.15 K). The Tm values of the other three parameters were detected at 338 K, a slightly higher temperature. These parameters had not improved in terms of the phase transition temperature. K

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 16. Snapshots of simulations using four parameters. The phases, Lc, Lβ, and Lα were predicted from the Tm values in the simulations. C, N, and O atoms in lipids are drawn as line model colored in gray, blue, and red, respectively. P atom is drawn as a vdW model colored in tan.

Nevertheless, the difference in the Tm value is ∼25 K, which was much smaller than those using GAFFlipid and Gaff. (Note that in the case using CHARMM, we used chamber software, which is in the AMBER software package, and not the original engine. If the original engine in the CHARMM software package is used, different results from this study may be obtained.) The simulation of DPPE with mod_Gaff2 and GROMACS indicated the same Tm value at 328 K, and the difference from the experimentally measured values was ∼10 K. The Tm values of lipid14 and CHARMM were 353 and 358 K, respectively. The differences from the experimentally measured values were slightly larger than those of mod_Gaff2 and GROMACS. Analyzing the AL values of DPPC, the profiles at the lowtemperature range differed according to the four parameters. The AL values of CHARMM and GROMACS were ∼50 Å2 at 243 K, whereas two AMBER-based parameters showed slightly smaller AL values, of which mod_Gaff2 and lipid14 were ∼45 and ∼42 Å2, respectively. At 273 K, the AL value of CHARMM gradually decreased, whereas those of the other three parameter sets gradually increased. Interestingly, after the main phase transition, the AL values were very similar in the four parameter sets. In the Lα phase, the alkyl chains completely melted and the fluidity of the lipids increased. According to our previous study, the vdW interaction in the alkyl chains strongly contributed to the thermal phase transition. In the Lα phase, the alkyl chains separated near the distance of the total vdW radius of carbon and hydrogen atoms by the vdW interaction of the alkyl chains. The total vdW radii of carbon and hydrogen (only carbon in GROMACS) in alkyl chains are very similar among four parameters. The AL values in the Lα phase, therefore, indicated similarity among the four parameters. Among the four parameters, the lipid conformations of DPPC in each phase differed slightly, as shown in Figure 16, and the conformations in the Lc and Lα phases were similar. However, in the Lβ phase, the lipid conformation of mod_Gaff2 differed from that of the other three parameters. In the previous

study, the alkyl-chain conformations showed a force in the opposite direction in the Lβ phase.13 The packing of alkyl chains changed from hexagonal to quasi-hexagonal shape. This conformational change occurred in the simulation using mod_Gaff2 but not with the other parameters, distinguishing mod_Gaff2 from the other parameters. It should be noted that to date we do not experimentally know the lipid conformation in the Lβ phase; therefore, we cannot validate the lipid conformation in the Lβ phase. However, according to the experimentally measured data, the phase transition appears between the Lc and the Lα phases. The conformation at this phase transition showed clearly that the alkyl chain packing changed from hexagonal to quasi-hexagonal shapes. We believe that the structure of the Lβ phase obtained by the simulation with mod_Gaff2 is similar to the original structure. The DPPE lipid conformation in the Lc phase with GROMACS is slightly different from those of the other parameters (Figure 16). In the temperature range from 273 to 313 K, the RF values of GROMACS were slightly larger than the other parameters. Therefore, most of the lipids, relative to the other parameters, had already melted before reaching the Tm. On the other hand, the lipid conformations in the Lα phase were similar among the four parameter sets. The alkyl chains completely melted, and the width of the bilayers shrank. These conformations did not differ significantly among the four parameters. From these results, the simulation for the lipids with the modified parameters largely reproduced the conformations with the other parameter sets. Moreover, the Tm values of DPPC and DPPE showed closer values to the experimentally measured ones than with the other parameters. To summarize, we have shown that the modified parameters can be improved for the description of the lipid thermal behaviors with a higher degree of agreement with experimentally accessible measures. L

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Simultaneous Analysis of Neutron and X-Ray Scattering Data. Biophys. J. 2008, 95, 2356−2367. (7) Schubert, T.; Schneck, E.; Tanaka, M. First Order Melting Transitions of Highly Ordered Dipalmitoyl Phosphatidylcholine Gel Phase Membranes in Molecular Dynamics Simulations with Atomistic Detail. J. Chem. Phys. 2011, 135, 055105. (8) Fidorra, M.; Heimburg, T.; Seeger, H. M. Melting of Individual Lipid Components in Binary Lipid Mixtures Studied by Ftir Spectroscopy, Dsc and Monte Carlo Simulations. Biochim. Biophys. Acta, Biomembr. 2009, 1788, 600−607. (9) Marrink, S. J.; de Vries, A. H.; Mark, A. E. Coarse Grained Model for Semiquantitative Lipid Simulations. J. Phys. Chem. B 2004, 108, 750−760. (10) Rodgers, J. M.; Sorensen, J.; de Meyer, F. J.; Schiott, B.; Smit, B. Understanding the Phase Behavior of Coarse-Grained Model Lipid Bilayers through Computational Calorimetry. J. Phys. Chem. B 2012, 116, 1551−1569. (11) Leekumjorn, S.; Sum, A. K. Molecular Studies of the Gel to Liquid-Crystalline Phase Transition for Fully Hydrated Dppc and Dppe Bilayers. Biochim. Biophys. Acta 2007, 1768, 354−365. (12) Heller, H.; Schaefer, M.; Schulten, K. Molecular Dynamics Simulation of a Bilayer of 200 Lipids in the Gel and in the Liquid Crystal Phase. J. Phys. Chem. 1993, 97, 8343−8360. (13) Ogata, K.; Uchida, W.; Nakamura, S. Understanding Thermal Phases in Atomic Detail by All-Atom Molecular-Dynamics Simulation of a Phospholipid Bilayer. J. Phys. Chem. B 2014, 118, 14353−14365. (14) Sonne, J.; Jensen, M. Ø.; Hansen, F. Y.; Hemmingsen, L.; Peters, G. H. Reparameterization of All-Atom Dipalmitoylphosphatidylcholine Lipid Parameters Enables Simulation of Fluid Bilayers at Zero Tension. Biophys. J. 2007, 92, 4157−4167. (15) Poger, D.; Van Gunsteren, W. F.; Mark, A. E. A New Force Field for Simulating Phosphatidylcholine Bilayers. J. Comput. Chem. 2010, 31, 1117−1125. (16) Klauda, J. B.; Brooks, B. R.; MacKerell, A. D.; Venable, R. M.; Pastor, R. W. An Ab Initio Study on the Torsional Surface of Alkanes and Its Effect on Molecular Simulations of Alkanes and a Dppc Bilayer. J. Phys. Chem. B 2005, 109, 5300−5311. (17) Klauda, J. B.; Monje, V.; Kim, T.; Im, W. Improving the Charmm Force Field for Polyunsaturated Fatty Acid Chains. J. Phys. Chem. B 2012, 116, 9424−9431. (18) Högberg, C.-J.; Nikitin, A. M.; Lyubartsev, A. P. Modification of the Charmm Force Field for Dmpc Lipid Bilayer. J. Comput. Chem. 2008, 29, 2359−2369. (19) Jorgensen, W. L.; Madura, J. D. Quantum and Statistical Mechanical Studies of Liquids. 25. Solvation and Conformation of Methanol in Water. J. Am. Chem. Soc. 1983, 105, 1407−1413. (20) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (21) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Mol. Graphics Modell. 2006, 25, 247−260. (22) 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, 9617. (23) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. (24) Case, D. A.; Cheatham, T. E., III; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M., Jr.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber Biomolecular Simulation Programs. J. Comput. Chem. 2005, 26, 1668−1688. (25) Koynova, R.; Caffrey, M. Phases and Phase Transitions of the Phosphatidylcholines. Biochim. Biophys. Acta 1998, 1376, 91−145. (26) 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, 10089−10092.

5. CONCLUSION We succeeded in improving the well-depth parameters to adapt with the Tm in DPPC. For the PC headgroups, simulations using the modified parameters reproduced Tm values similar to experimentally measured ones. The errors of the Tm values were around 15 K. These error values were still well improved when compared with those from the simulation using GAFFlipid and Gaff. The lipids in the simulations using the modified parameters were able to keep their structural properties like the simulation using the original parameters. The simulation using mod_Gaff2 showed similar results as other parameter sets, lipid14, CHARMM, and GROMACS. We therefore conclude that our improved parameter is very efficient for simulation, not only for the lipid bilayers but also for the transmembrane proteins embedded in the transmembrane.



ASSOCIATED CONTENT

S Supporting Information *

Atom names of lipids, detailed graphs of AL, RF, SCD, lateral radial distribution function, electron density profile, and tables of atomic charges in the lipids. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b01656.



AUTHOR INFORMATION

Corresponding Author

*Phone: +81-(0)48467-9477. Fax: +81-(0)48467-8503. E-mail: [email protected]. Author Contributions

K.O. performed MD simulations, analyzed the results, wrote the paper, and was involved in the design of the study. S.N. designed the study, handled overall project strategy and management, and wrote the paper. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We thank the KAITEKI Institute (Mitsubishi Chemical Holdings) for their support and Dr. T. Iitaka of RIKEN and Ms W. Uchida of Tokyo Institute of Technology for fruitful discussions. Some of the calculations in this study were performed using the RIKEN Integrated Cluster of Clusters (RICC) Facility.



REFERENCES

(1) Voet, J. G.; Voet, D. Biochemistry; John Wiley& Sons: Hoboken, NJ, 2004. (2) Chen, S. C.; Sturtevant, J. M.; Gaffney, B. J. Scanning Calorimetric Evidence for a Third Phase Transition in Phosphatidylcholine Bilayers. Proc. Natl. Acad. Sci., U.S.A. 1980, 77, 5060−5063. (3) Lewis, R. N.; Mak, N.; McElhaney, R. N. A Differential Scanning Calorimetric Study of the Thermotropic Phase Behavior of Model Membranes Composed of Phosphatidylcholines Containing Linear Saturated Fatty Acyl Chains. Biochemistry 1987, 26, 6118−6126. (4) Nagle, J. F.; Tristram-Nagle, S. Structure of Lipid Bilayers. Biochim. Biophys. Acta 2000, 1469, 159−195. (5) Kusube, M.; Matsuki, H.; Kaneshina, S. Thermotropic and Barotropic Phase Transitions of N-Methylated Dipalmitoylphosphatidylethanolamine Bilayers. Biochim. Biophys. Acta 2005, 1668, 25−32. (6) Kučerka, N.; Nagle, J. F.; Sachs, J. N.; Feller, S. E.; Pencer, J.; Jackson, A.; Katsaras, J. Lipid Bilayer Structure Determined by the M

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (27) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (28) Loncharich, R. J.; Brooks, B. R.; Pastor, R. W. Langevin Dynamics of Peptides: The Frictional Dependence of Isomerization Rates of N-Acetylalanyl-N′-Methylamide. Biopolymers 1992, 32, 523− 535. (29) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (30) Dickson, C. J.; Madej, B. D.; Skjevik, Å. A.; Betz, R. M.; Teigen, K.; Gould, I. R.; Walker, R. C. Lipid14: The Amber Lipid Force Field. J. Chem. Theory Comput. 2014, 10, 865−879. (31) Brooks, B. R.; Brooks, C. L., III; Mackerell, A. D., Jr.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. Charmm: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545−1614. (32) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. Gromacs 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (33) Crowley, M. F.; Williamson, M. J.; Walker, R. C. Chamber: Comprehensive Support for Charmm Force Fields within the Amber Software. Int. J. Quantum Chem. 2009, 109, 3767−3772. (34) Poger, D.; Mark, A. E. On the Validation of Molecular Dynamics Simulations of Saturated and Cis-Monounsaturated Phosphatidylcholine Lipid Bilayers: A Comparison with Experiment. J. Chem. Theory Comput. 2009, 6, 325−336. (35) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Intermolecular Forces; Reidel: Dordrecht, 1981. (36) Parrinello, M.; Rahman, A. Crystal Structure and Pair Potentials: A Molecular-Dynamics Study. Phys. Rev. Lett. 1980, 45, 1196−1199. (37) Seelig, J.; Niederberger, W. Deuterium-Labeled Lipids as Structural Probes in Liquid Crystalline Bilayers. Deuterium Magnetic Resonance Study. J. Am. Chem. Soc. 1974, 96, 2069−2072. (38) Seelig, J.; Waespe-Sarcevic, N. Molecular Order in Cis and Trans Unsaturated Phospholipid Bilayers. Biochemistry 1978, 17, 3310−3315. (39) Petrache, H. I.; Dodd, S. W.; Brown, M. F. Area Per Lipid and Acyl Length Distributions in Fluid Phosphatidylcholines Determined by 2h Nmr Spectroscopy. Biophys. J. 2000, 79, 3172−3192. (40) Kučerka, N.; Nieh, M.-P.; Katsaras, J. Fluid Phase Lipid Areas and Bilayer Thicknesses of Commonly Used Phosphatidylcholines as a Function of Temperature. Biochim. Biophys. Acta: Biomembr. 2011, 1808, 2761−2771. (41) Kučerka, N.; Tristram-Nagle, S.; Nagle, J. Structure of Fully Hydrated Fluid Phase Lipid Bilayers with Monounsaturated Chains. J. Membr. Biol. 2006, 208, 193−202. (42) Rappolt, M.; Hickel, A.; Bringezu, F.; Lohner, K. Mechanism of the Lamellar/Inverse Hexagonal Phase Transition Examined by High Resolution X-Ray Diffraction. Biophys. J. 2003, 84, 3111−3122. (43) Leftin, A.; Brown, M. F. An Nmr Database for Simulations of Membrane Dynamics. Biochim. Biophys. Acta: Biomembr. 2011, 1808, 818−839. (44) Douliez, J. P.; Leonard, A.; Dufourc, E. J. Restatement of Order Parameters in Biomembranes: Calculation of C-C Bond Order Parameters from C-D Quadrupolar Splittings. Biophys. J. 1995, 68, 1727−1739. (45) Leftin, A.; Molugu; Trivikram, R.; Job, C.; Beyer, K.; Brown; Michael, F. Area Per Lipid and Cholesterol Interactions in Membranes from Separated Local-Field 13c Nmr Spectroscopy. Biophys. J. 2014, 107, 2274−2286. (46) Nagle, J. F.; Zhang, R.; Tristram-Nagle, S.; Sun, W.; Petrache, H. I.; Suter, R. M. X-Ray Structure Determination of Fully Hydrated L Alpha Phase Dipalmitoylphosphatidylcholine Bilayers. Biophys. J. 1996, 70, 1419−1431. (47) Tristram-Nagle, S.; Petrache, H. I.; Nagle, J. F. Structure and Interactions of Fully Hydrated Dioleoylphosphatidylcholine Bilayers. Biophys. J. 1998, 75, 917−925.

N

DOI: 10.1021/acs.jpcb.5b01656 J. Phys. Chem. B XXXX, XXX, XXX−XXX