Biomolecular Simulations with the Transferable Potentials for Phase

Article pubs.acs.org/JPCB

Biomolecular Simulations with the Transferable Potentials for Phase Equilibria: Extension to Phospholipids Navendu Bhatnagar,† Ganesh Kamath,‡ and Jeffrey J. Potoff*,† †

Department of Chemical Engineering and Materials Science, Wayne State University, Detroit, Michigan 48202, United States Department of Chemistry, University of Missouri−Columbia, Columbia, Missouri 65211-7600, United States

‡

S Supporting Information *

ABSTRACT: The Transferable Potentials for Phase Equilibria (TraPPE) is extended to zwitterionic and charged lipids including phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylserine (PS), and phosphatidylglycerol (PG). The performance of the force field is validated through isothermal−isobaric ensemble (NPT) molecular dynamics simulations of hydrated lipid bilayers performed with the aforementioned head groups combined with saturated and unsaturated alkyl tails containing 12−18 carbon atoms. The effects of water model and sodium ion parameters on the performance of the lipid force field are determined. The predictions of the TraPPE force field for the area per lipid, bilayer thickness, and volume per lipid are within 1−5% of experimental values. Key structural properties of the bilayer, such as order parameter splitting in the sn-2 chain and X-ray form factors, are found to be in close agreement with experimental data.

1. INTRODUCTION Traditional approaches to the parametrization of nonbonded interactions for biomolecular force fields involve the optimization of parameters to reproduce condensed phase properties, such as liquid densities and heats of vaporization at ambient conditions.1−4 Additional input from quantum mechanical calculations, such as binding energies with water3,5 or rare gas atoms,6,7 may be used to further optimize partial charge distributions and Lennard-Jones parameters. Free energies of solvation in water and other solvents, such as hexane, cyclohexane, or 1-octanol, have also been used as a target parameters in the optimization process.4,9,10 On the other end of the spectrum are force fields for the prediction of vapor−liquid equilibria.11−15 These force fields typically make use of quantum chemical calculations for the determination of partial charge distributions, while LennardJones parameters are optimized to reproduce complete vapor liquid coexistence curves from the normal boiling point to the critical point. Additional data such as vapor pressures16,17 or mixture phase equilibria may also be used in the optimization process.18,19 By optimizing nonbonded parameters over a broad range of temperatures and pressures, these force fields are capable of making reliable predictions for fluid physical properties over a broad range of temperatures, pressures, and phases (gas, liquid, and solid), for pure fluids and complex mixtures. In recent years, generalized versions of biomolecular force fields have been published,21,22 and these have been used to provide valuable insight for a broad array of chemical in addition to biological systems. On the other hand, there has been little crossover for force fields developed for the prediction of vapor−liquid equilibria into the realm of © 2013 American Chemical Society

biomolecular computation. In the past, this was due to limitations of the parameter sets; however, in recent years, parameter sets have expanded to the point where the simulation of biological systems is now a possibility. Lipid bilayers composed of phospholipids are crucial to sustaining life, forming the boundaries of cells and acting as scaffolds to support a variety of proteins. Because of their central importance in understanding life processes, lipid bilayers have been studied extensively, through both experiments and computer simulations. Early attempts to simulate lipid bilayers relied on “ad-hoc” force fields, with parameters derived from analogous compounds and combined to form the lipids of interest. An example of this is the well-known “Berger” force field for phosphatidylcholine, which was derived by combining bonded interactions from GROMOS with nonbonded interactions from the Optimized Potentials for Liquid Simulations United-Atom (OPLS-UA) force field23 and partial charges from the work of Chiu et al.24 This model included new parameters for the alkyl tail, which were optimized to reproduce the heat of vaporization of pentadecane. Additional unitedatom force fields for phoshatidylcholine lipids have appeared in the literature based on the GROMOS force field, of which there are numerous variants. Chiu et al. combined partial charge distributions from a prior work,24 with Lennard-Jones parameters optimized by fitting to reproduce liquid densities and heats of vaporization of model compounds to produce the 43A1-S3 parameter set for lipids.25 While providing a good reproduction of bilayer properties, such as area per lipid, the Received: May 1, 2013 Revised: July 24, 2013 Published: July 29, 2013 9910

dx.doi.org/10.1021/jp404314k | J. Phys. Chem. B 2013, 117, 9910−9921

The Journal of Physical Chemistry B

Article

also provides the foundation for another lipid force field known as Stockholm Lipids (Slipids).44,45 Given the extensive data that exist in the literature, phospholipid bilayers are the ideal system to assess the reliability of intermolecular potentials developed for vapor− liquid equilibria calculations when used in biological systems. In this work, the Transferable Potentials for Phase Equilibria (TraPPE)11,46−53 are extended to phospholipids that contain phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylserine (PS), and phosphatidylglycerol (PG) head groups and saturated and unsaturated alkyl tails with 12−18 carbon atoms. To highlight the transferable nature of the TraPPE force field, Lennard-Jones parameters developed for the prediction of vapor−liquid equilibria for low molecular weight organic molecules are used without modification. Partial charges are determined from a CHELPG analysis of HF/631+G(d,p) ab initio calculations, while torsional parameters are refit to match the rotational barriers predicted by the CHARMM C36 lipid force field and MP2/6-31G(d) ab initio calculations. A variety of quantities, including the area per lipid, volume per lipid, bilayer thickness, electron density, and X-ray form factors are calculated from molecular dynamics simulations in the tensionless NPT ensemble and are shown to be in close agreement with experimental data.

43A1-S3 force field introduced subtle inconsistencies in the nonbonded parameters used for lipids and other compounds. A self-consistent GROMOS lipid force field for DPPC based on the G53A6 parameter set, referred to as G53A6L, appeared shortly thereafter.26 Additional models have been proposed for phosphatidylglycerol (PG) lipids by Dickey and Faller who combined nonbonded parameters from the Berger lipid force field with GROMOS8727 and Kukol,28 who combined the G53A6 Lennard-Jones parameters with partial charges from Zhao et al.29 A revised set of bonded parameters has been proposed for the Kukol POPG force field that improve its performance with respect to the prediction of deuterium order parameters.30 Development of lipid force fields based on the AMBER parameter set followed a path similar to GROMOS. Parameters from the Generalized Amber Force Field (GAFF) were combined with RESP derived partial charges31 to produce a model for POPC. However, simulations in the tensionless NPT ensemble led to bilayer condensation, with area per lipid values approximately 20% below experiment. For this parameter set, an applied surface tension of 30 mN/m was required to reproduce the experimental area per lipid.32 LIPID11 uses a modular approach, analogous to proteins, to build lipid force fields from GAFF, although the resulting models require simulation using constant area, or with a nonzero surface tension to prevent condensation above the gel phase transition temperature.33 The tendency for bilayer condensation was eliminated by refitting of the Lennard-Jones parameters for carbon and hydrogen atoms in the alkyl tails to reproduce the heat of vaporization and liquid density for pentadecane.34 This parametrization is known as GAFFlipid. Like GAFF, CHARMM is an all-atom force field for the simulation of nucleic acids, proteins, carbohydrates, and lipids. The CHARMM lipid force field has a long history and has undergone continuous revision as computer hardware advances have provided access to the simulation of larger systems and longer time scales. While lipid bilayer simulations were performed with the C22 parameters set for approximately 1 ns,35,36 the latest C36 parameter set was validated with simulations of over 100 ns.37 For simulations of DPPC with the C22 parameters, NPT simulations of approximately 1 ns showed good agreement with experimental area per lipid values. However, later simulations of DPPC bilayers in the NPγT ensemble revealed that for γ = 0 the C22 parameter set produced area per lipid values that were approximately 14% lower that experimental values.38 The C27 parametrization involved the optimization of Lennard-Jones parameters for the alkyl tails, and torsional potentials for the phosphate headgroup.39 Despite these improvements, the C27 parameter set still exhibited a tendency toward condensation at temperatures above the gel phase transition temperature when used in tensionless NPT simulations.40 Through reoptimization of the partial charges for the phospholipid headgroup using the RESP method,41,42 Sonne et al. produced a variant of the C27 force field that was capable of predicting an area per lipid for DPPC within 6% of experimental values in NPT simulations.43 The C36 parameter set included significant reparameterization of torsional potentials, partial charge distributions, and LennardJones parameters for both the lipid headgroup and the alkyl tails. Extensive NPT simulations have shown that this latest parametrization reproduces closely a broad range of experimental data, including area per lipid.37 The C36 parameter set

2. FORCE FIELD Nonbonded interactions in the TraPPE force field are described by pairwise-additive 12-6 Lennard-Jones potentials and Coulombic interactions of partial charges ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ qiqj σij σij Uij = 4εij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ + r ⎝ rij ⎠ ⎦ 4πε0rij ⎣⎝ ij ⎠

(1)

where rij, εij, σij, qi, and qj are the separation, LJ well depth, pseudoatom diameter, and partial charges, respectively, for the pair of interaction sites i and j and ε0 is the permittivity of vacuum. A united-atom representation is used for all CHx groups; i.e., hydrogen atoms bonded to carbon atoms are not represented explicitly and are instead combined with the carbon atoms to form a single interaction action site or “pseudo-atom”. Lennard-Jones parameters for each interaction site were taken from previous TraPPE publications and used without modification.11,46,48,49,51,52,54 Parameters for interaction sites of different types were determined with the standard Lorentz− Berthelot combining rules. Partial charges were determined from ab initio calculations performed on complete head groups, but without the hydrocarbon tails. Previous calculations have shown splitting of the headgroup into fragments may result in partial charge distributions that lead to area per lipid values that are significantly lower than experimental values.55 Consistent with the TraPPE force field, CHX groups in the hydrocarbon tail were given a partial charge of 0. For zwitterionic lipids (PC, PE), structures were optimized using Hartree−Fock theory with the 6-31G(d,p) basis set, while for charged lipids (PS, PG) the 6-31+G(d,p) basis set was used. Partial charges were extracted from the electrostatic potential energy surfaces using the CHELPG (CHarges from ELectrostatic Potentials using a Grid based method) methodology.56 Charges on symmetrical or identical interaction sites were averaged to simplify the model. All quantum calculations were performed with Gaussian 09.57 9911



Article

Figure 1. Schematic of model compounds used in this study for optimization of dihedral parameters: (a) choline, (b) dimethyl phosphate (DMP), (c) 2-hexene, (d) hexane, (e) esterified glycerolphosphate analogue ((M)PGLY), and (f) esterified glycerol analogue ((M)EGLY).

Each phospholipid was modeled as a fully flexible molecule, with conformational behavior controlled by a combination of harmonic potentials for bond stretching and angle bending

U = k(x − x0)2

dihedral potential for the n-alkane tail with the TraPPE nonbonded parameters resulted in a significant overprediction of trans relative to cis conformations. Resolving this issue required the optimization of new parameters for the dihedral potential that took into account explicit 1-4 Lennard-Jones and Coulombic interactions. Lipid molecules were split into a series of model compounds, which are shown in Figure 1. Choline and dimethyl phosphate were used to parametrize the torsions in lipid headgroup, ethyl acetate for glycerol linker, and hexane and hexene for the tail region. Dihedral parameters derived from these compounds were further optimized in calculations of rotational barriers in (M)EGLY and (M)PGLY. The use of larger model compounds for the optimization of dihedral parameters was found to be crucial to reproducing experimental order parameter data, including the phenomena of order parameter splitting near the headgroup. The branched region in (M)EGLY formed by atoms C1−C2−C3−O31 and C1−C2−O21−C21 contains multiple redundant dihedrals, the parameters for which were fit simultaneously. For each dihedral of interest, rotational barriers were determined from relaxed potential energy scans performed with the CHARMM software, version 34. For the lipids containing the phosphatidylcholine headgroup, parameters were optimized to reproduce rotational barriers given by the CHARMM 36 force field.37 The CHARMM C36 force field was selected as the target data over quantum calculations since these parameters have already been optimized to reproduce QM data, and represent the current “state of the art” for allatom lipid simulations. Additional rotational barriers found in phosphatidylserine and phosphatidylglycerol were optimized to reproduce data extracted from relaxed potential energy scans from MP2/6-31G(d) quantum chemical calculations. In Figure 2, a comparison between selected rotational barriers predicted with the optimized TraPPE parameters and the CHARMM C36 force field is presented. Data for rotational

(2)

where U corresponds to the bond stretching or bending energy and x and x0 are the instantaneous (bond length or angle) and equilibrium (bond length or angle), respectively. Traditionally, the TraPPE force field uses fixed bond lengths; however, in this work, flexible bonds were used to simplify the implementation of the force field in molecular dynamics simulations. For some bond angles, an additional Urey−Bradley term was included. Values of the bond stretching and bending constants were taken from the CHARMM C36 force field.37 Rotations around dihedral angles are described with a cosine series m

Utors =

∑ ki[1 + cos(iφ − δ)] i=1

(3)

where ki are force constants, n is the periodicity or multiplicity, φ is the dihedral angle, and δ the phase shift. In prior implementations of the TraPPE force field, 1-4 Lennard-Jones and electrostatic interactions were not calculated explicitly. Instead, the cosine series was used to represent the total energy for barriers to rotation. However, this poses problems for the use of TraPPE in most molecular dynamics simulation engines, which follow the conventions of biomolecular force fields, such as CHARMM and AMBER, and explicitly include 1-4 Lennard-Jones and electrostatic interactions as part of the calculation of rotational barriers. Naively combining the CHARMM torsional potentials with the TraPPE Lennard-Jones and partial charges produced rotational barriers with significant errors compared to quantum chemical calculations. For example, the combining of the CHARMM 9912



Article

interactions of the bilayer with itself through periodic boundary conditions. 128 Na+ were added to the PS and PG systems to maintain charge neutrality. This corresponds to a salt concentration of approximately 60 mM. Interactions between water molecules were governed by the SPC/E force field.58 Simulations were performed with NAMD 2.9 in the isothermal−isobaric ensemble.59 Lennard-Jones interactions were truncated at 14 Å, with a smooth switching function starting at 12 Å. It should be noted that the TraPPE force field was originally developed using an analytical tail correction, which was neglected here for simple implementation in typical molecular dynamics simulation engines. The choice of a 14 Å cutoff was based on prior calculations which have shown truncation of the Lennard-Jones potential at 14 Å produces liquid densities that are in within 1.5% of the original tailcorrected potential.60 Electrostatic interactions were calculated using the particle mesh Ewald (PME). Temperature was controlled using Langevin dynamics with a damping coefficient of 5.0 ps−1. The pressure was maintained at 1.013 bar with the Langevin piston Nosé−Hoover method,61,62 with an oscillation period of 100 fs and piston decay period of 50 fs. The cell volume was allowed to fluctuate independently in the x−y plane and the z-coordinate axis. Initial coordinates for each system were generated using CHARMM-GUI (http://www. charmm-gui.org/). All the systems were initially minimized to the lowest energy state, followed by slow heating to a temperature of 323 K for DPPC, and 303 K for all other lipids. These temperatures were chosen to be above the gel phase transition temperatures of the respective lipids and consistent with available experimental data. Systems were simulated for 100 ns, using a time step of 2.0 fs.

Figure 2. Rotational barriers predicted for various dihedrals in phosphatidylcholine-based lipids using the TraPPE (red circles) and CHARMM C36 force fields (black line): (a) N−C12−C11−O12 (choline), (b) C1−C2−C3−O31 (M)EGLY, (c) P−O11−C1−C2 (M)PGLY, (d) CT2−CT2−CT2−CT2 (hexane), (e) O12−P−O11− C1 (M)PGLY, (f) O21−C21−C22−C23 (M)EGLY, (g) C29 C210−C211−C212 (2-hexene), and (h) O13−P−O11−C1 (M)PGLY.

barriers not presented here may be found in the Supporting Information (Figures S1, S2, and S3). Overall, the data show close agreement with predictions of the C36 parameter set. Subtle differences were observed for some of the rotational barriers in (M)EGLY and (M)PGLY. For the O21−C21− C22−C23 dihedral in (M)EGLY, the TraPPE force field overpredicts the magnitude of the barrier at 110° by approximately 0.6 kcal/mol, while for the C1−C2−C3−O31 dihedral, the locations of minima are shifted by 2−3° compared to the C36 parameter set. These small differences arise due to limitations in the united-atom representation, which does not provide as much fine-grained control over 1-4 Lennard-Jones interactions as an all-atom representation.

4. RESULTS AND DISCUSSION 4.1. Area per Lipid. The area per lipid is a key structural property of lipid bilayers and is governed by a complex interplay of interactions between lipids and water. Of particular importance are interactions between water and oxygen atoms in the glycerol linker region. If these interactions are too weak, bilayers exhibit a tendency in tensionless NPT simulations to condense to area per lipid values indicative of the gel phase. This was illustrated clearly in simulations of DPPC bilayers using the GROMOS 53A6 parameter set with two sets of charges, GROMOS96, and those taken from Chiu et al.24,26 Simulations using the GROMOS96 charge set produced an area per lipid of 42.6 Å2, while those using the Chiu charges produced an area per lipid of 63.1 Å2. The CHARMM C27 force field for DPPC also displayed a propensity to condense at temperatures above the gel-phase transition temperature.43 Application of new partial charges determined from a RESP fit of RHF/6-31G(d) ab initio data to the C27 parameter set produced a model capable of reproducing the experimental area per lipid.43 For both the GROMOS and CHARMM C27 force fields, the revised partial charge distributions feature significant increases in the magnitude of the partial charges for oxygen atoms in the fatty acid section of the lipid tails.26,43,63 Because of the importance of water interactions with the lipid headgroup, the choice of water model can have a significant impact on the predicted values of the area per lipid. Care must be taken to ensure the derived partial charge distribution for the lipids of interest is compatible with the chosen water model. In this work, partial charges were derived from a CHELPG analysis of HF/6-31G(d,p) and HF/6-31+G(d,p) ab initio calculations, which has been found by our group to produce

3. COMPUTATIONAL DETAILS The TraPPE force field was developed for four lipid head groups, including phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylserine (PS), and phosphatidylglycerol (PG), combined with saturated and unsaturated alkyl tails containing 12−18 carbon atoms. The performance of the force field was validated with simulations performed on 10 lipids: dipalmitoylphosphatidylcholine (DPPC), dimyristoylphosphatidylcholine (DMPC), dilauroylphosphatidylcholine (DLPC), dioleoylphosphatidylcholine (DOPC), 1-palmitoyl-2-oleoylphosphatidylcholine (POPC), 1-palmitoyl-2-oleoyl-phosphatidylethanolamine (POPE), dioleoylphosphatidylserine (DOPS), dilauroylphosphatidylglycerol (DLPG), dimyristoylphosphatidylglycerol (DMPG), and dioleoylphosphatidylglycerol (DOPG). It is expected that the force fields described here will be transferable to other combinations of head and tail groups not explicitly defined in this work. Calculations were performed on bilayers containing 128 lipids, solvated by approximately 50 waters per lipid for zwitterionic lipids, and 90 waters per lipid for charged lipids. Additional waters were included for charged lipids to provide sufficient hydration of ions used for charge neutralization, and minimize the 9913



Article

performed for 100 ns with the Roux Na+ parameters. The final configurations from these simulations were used as input to simulations performed with the Aqvist parameters for Na+. The evolution of the area per lipid during the initial 100 ns simulation with the Roux Na+ parameters and after the change to the Aqvist Na+ parameters is shown in Figure 4. Equilibration of the area per lipid occurred within 20 ns of the introduction of the Aqvist parameters for Na+. The predicted area per lipid values are shown in Figure 3 versus experimental data and the numerical values are listed in Table 1. Simulations performed using the Roux Na+ ion produced area per lipid values for PG containing lipids that were approximately 3 Å2 lower than experimental values. This was due to extensive binding of Na+ to glycerol and phosphate oxygen atoms in neighboring lipid head groups, as shown in Figure 5. Cation binding to the lipid head groups reduces the hydration of the bilayer slightly, leading to bilayer contraction. The Aqvist parameters for Na+ produce a minimum in the interaction energy for Na+−oxygen interactions that is about 10 kcal/mol weaker than the Roux parameters, which is shown in Figure S7 in the Supporting Information. Comparison of radial distribution functions, Figure 6, illustrates how the difference in interaction strength lead to weaker association of the Aqvist Na+ with various oxygen atoms compared to the Roux Na+ parameters. While the location of the first peak was the same for both parameter sets, the peak height was reduced substantially for the Aqvist Na+ model. In Table 2, the fraction of Na+ bound to the various oxygen atoms in DLPG is listed. Na+ within 3 Å of an oxygen atom were considered bound; this distance was determined from the location of the first minimum in the Na+−oxygen radial distribution functions. The data show that, for the Roux parameters, approximately 59% of Na+ was bound to the lipid headgroup. For the Aqvist parameters, the fraction of bound Na+ dropped to approximately 45%. All of the Na+ binding occurred at the sp2-hybridized oxygens in the phosphate and carboxyl groups. The reduction in ion binding leads to greater hydration of the phospholipid head groups. This is illustrated in the electron densities calculated for DLPG, Figure 7, which show that greater hydration of the bilayer is achieved with the Aqvist Na+ parameters. Greater hydration causes an expansion of the bilayer surface area by approximately 3 Å2, leading to improved predictions of the area per lipid. Simulations using the Aqvist parameters for Na+ produced area per lipid values with an AAD of 0.6 Å2, compared to an AAD of 3.1 Å2 observed in simulations using the Roux Na+ parameters. This compares favorably to the predictions the Slipids all-atom force field for PG-containing lipids (AAD = 1.07 Å2), which also used the Aqvist parameters for Na+.68 From fluctuations in the area per lipid, is possible to estimate the area compressibility modulus KA

accurate free energies of hydration for solutes modeled by the TraPPE force field in SPC/E water.60 Predictions of the TraPPE force field for area per lipid are plotted versus experiment in Figure 3, and numerical data are listed in Table 1.

Figure 3. Area per lipid predicted by the TraPPE force field for lipid bilayers composed of phosphatidylcholine (red circles), phosphatidylethanolamine (orange squares), phosphatidylserine (blue diamonds), and phosphatidylglycerol (green triangles).

The evolution of the area per lipid with time is presented in Figure 4 for each lipid. For lipids that contain the phosphatidylcholine headgroup (DPPC, DMPC, DLPC, DOPC, POPC), the TraPPE force field predicted area per lipid values with an absolute average deviation of 1.02 Å2, which is of similar accuracy as the most accurate all-atom force fields, CHARMM C36,37 and Slipids.44 For POPE, TraPPE predicts an area per lipid of 62.3 Å2, which compares favorably to the experimental value of 59 Å2.64 To understand the effect of water model on the predicted area per lipid, additional calculations were performed for DPPC using the TIP3P force field. The use of TIP3P produced a slight increase in the predicted area per lipid to 68.3 Å2, compared to 63.6 Å2 for simulations using SPC/E. Density profiles show slightly deeper penetration of TIP3P water into the lipid bilayer compared to simulation performed with SPC/E water (Figure S5 in the Supporting Information). Radial distribution functions for water interacting with oxygens in the glycerol region (Figure S6 in the Supporting Information) show a measurable increase in the first peak for TIP3P water compared to SPC/E, illustrating how stronger water−lipid interactions lead to an increase in hydration of the phospholipid headgroup and greater area per lipid. Simulation of systems containing charged lipids, such as those containing phosphatidylserine or phosphatidylglycerol, is even more complex, since they require a neutralizing cation, such as Na+ or K+. Numerous simulations have shown that cations may bind to phospholipid head groups of zwitterionic lipids, leading to significant reductions in the predicted area per lipid. For example, simulations of POPC, DPPC showed that the addition of Na+ or K+ caused a reduction in the area per lipid by 10−15%.65−67 Therefore, as in the case of water, it is important that the parametrization of the cations be consistent with the phospholipid model. To understand the effect of cation parameters on the predicted area per lipid, calculations were performed on phosphatidylserine and phosphatidylglycerol containing lipids with 128 Na+ as the counterion using parameters from Roux8 or Aqvist20 for the sodium ion. NPT simulations were initially

KA =

2AL βnLσ 2

(4)

where β = 1/kbT, nL is the number of lipids, and σ2 is the variance in the area per lipid AL. These data are presented in Table 1. For PC-containing lipids, TraPPE KA values are approximately 60% larger than experiment. A similar overprediction was observed in simulations of the GROMOS 53A6L force field.69 This could be due to system size effects, and periodic boundary conditions, which suppress membrane undulations.70,71 Recent calculations performed with GAFFlipid 9914



Article

Table 1. Area per Lipid, Volume per Lipid, Bilayer Thickness, and Area Compressibility Modulus Predicted by TraPPE United Atom Force Fielda Alipid (Å2) expt74,86 TraPPE CHARMM3637 GAFFlipid34 Slipids44 GROMOS 53A6L26,87 Berger23,25 43A1-S325

63.1 63.6 ± 0.6 62.9 61.2 62.6 63.1 65.5 63.7

expt82,88 TraPPE CHARMM3637 GAFFlipid34 Slipids44 GROMOS 53A669 Berger23,25 43A1-S325

60.6 61.4 ± 0.5 60.8 59.9 60.8 61.6 62.6 62.1

expt82 TraPPE CHARMM3637 GAFFlipid34 Slipids44 GROMOS 53A669 Berger23,25 43A1-S325

63.2 62.8 ± 0.4 64.4 59.9 62.4 63.2 65.4 63.0

expt74,89 TraPPE CHARMM3637 GAFFlipid34 Slipids44 GROMOS 53A669 Berger23,25 43A1-S325

67.4 68.8 ± 0.5 69.0 66.8 68.0 64.9 66.3 66.0

expt90,91 TraPPE CHARMM3637 GAFFlipid34 Slipids44 GROMOS 53A669

68.3 66.3 ± 0.6 64.7 63.7 64.6 63.8

expt64,92 TraPPE CHARMM3637 GAFFlipid34 Slipids44

59 62.3 ± 0.4 59.2 55.64 56.3

expt81 TraPPE/Na+ (Roux) TraPPE/Na+ (Aqvist) Slipids68

65.3 67.7 ± 0.5 70.7 ± 0.4 64.1


65.6 61.5 ± 0.7 65.5 ± 0.1 64.2

Vlipid (Å3) DPPC (T = 323 K) 1229 1231 ± 0.8 − 1265.4 1201 1226 1226 1209 DMPC (T = 303 K) 1101 1106 ± 0.6 − 1117.8 1060 1077 1112 1086 DLPC (T = 303 K) 991 996 ± 0.9 − 1008.8 951 969 1013 977 DOPC (T = 303 K) 1303 1310 ± 0.7 − 1327.4 1262 1284 1343 1291 POPC (T = 303 K) 1256 1264 ± 1.0 − 1277 1213 1232 POPE (T = 303 K) 1175 1228 ± 0.9 − 1185.21 1153 DOPS (T = 303 K) 1228 1251 ± 1.1 1260.9 ± 7.5 1222 DLPG (T = 303 K) 953.6 938.4 ± 1.2 946.4 ± 4.3 922

9915

DHH (Å)

DB (Å)

KA (mN/m)

38.0 37.5 ± 0.2 − 37.6 37.7 35.7 34.7 35.7

39.0 38.7 ± 0.3 − − 38.3 38.3 37.3 38.0

231 371 207 274 238 414 − −

± ± ± ± ± ±

35.3 34.5 ± 0.1 − 33.6 34.5 32.7 32.3 32.6

36.3 36.0 ± 0.2 − − 35.3 34.9 35.5 35.0

234 383 − 299 250 475 − −

± 23 ± 20

30.8 30.5 ± 0.3 − 31.6 30.1 28.5 27.8 28.5

31.4 31.7 ± 0.2 − − 30.4 30.7 31.0 31.0

− 468 − 291 268 461 − −

36.7 37.1 ± 0.1 − 37.6 36.6 36.3 37.2 36.6

38.7 38.1 ± 0.3 − − 37.2 38.9 40.4 39.1

254 406 − 314 256 389 − −

37 37.1 ± 0.3 − 37.6 36.5 34.6

36.8 38.1 ± 0.2 − − 38.5 38.7

180−330 300 ± 16 − 391 ± 81 298 ± 30 404 ± 55

40.0 37.0 ± 0.1 − 43.4 41.1

− 39.5 ± 0.3 − − 41.6

233 264 ± 24 − 484 ± 34 282 ± 29

38.4 38.0 ± 0.2 36.9 ± 0.3 38.1

37.6 37.8 ± 0.2 36.1 ± 0.6 38.5

− 378 ± 56 311 ± 30 261 ± 27

− 30.2 ± 0.2 29.0 ± 0.3 27.0

29.1 31.1 ± 0.2 29.2 ± 0.1 29.0

− 121 ± 20 311 ± 23 165 ± 30

20 40 14 22 35 59

± 75 ± 29 ± 10

± 86 ± 64 ± 24 ± 96

± 27 ± 39 ± 29 ± 19



Article

Table 1. continued Alipid (Å2)

a

expt93 TraPPE/Na+ (Roux) TraPPE/Na+ (Aqvist) Slipids68 CHARMM94

65.1 61.3 ± 0.5 65.1 ± 0.1 63.3 58.5


70.8 68.4 ± 0.2 72.2 ± 0.4 70.8

Vlipid (Å3) DMPG (T = 303 K) 1057.4 1045.9 ± 1.3 1054 ± 8.2 1030 − DOPG (T = 303 K) 1265 1252.5 ± 1.1 1260.1 ± 9.7 1232

DHH (Å)

DB (Å)

KA (mN/m)

− 34.1 ± 0.1 32.3 ± 0.2 31.8 −

32.5 34.4 ± 0.2 32.4 ± 0.1 32.4 −

− 103 ± 15 242 ± 18 161 ± 35 −

− 36.8 ± 0.2 35.4 ± 0.1 34.5

35.7 36.6 ± 0.1 34.9 ± 0.4 35.0

− 200 ± 37 278 ± 32 241 ± 36

Alipid, area per lipid; Vlipid, volume per lipid; DHH, bilayer thickness; DB, Luzzati thickness; KA, area compressibility modulus.

Figure 4. Evolution of the area per lipid (black line) and the average area per lipid (green dashed line) predicted by the TraPPE force field. Experimental data (red line). For PG-containing lipids, simulations performed with the Roux Na+ parameters: evolution of the area per lipid (black line) and average area per lipid (green-dashed line); simulations performed with the Aqvist Na+ parameters: evolution of the area per lipid (orange line) and the average area per lipid (blue dashed line).

for a variety of system sizes suggest KA is sensitive to system size, and that simulations of significantly larger systems may be necessary to predict KA accurately from fluctuations in the area per lipid observed during simulation.34,71 The predicted area compressibility modulus for PG lipids is significantly lower than for PC lipids, and follows the same trend as the Slipids force field.68 This is the result of hydrogen bonding between lipids, as well as binding of Na+ ions to the lipid head groups, which lead to a significantly larger variance in the area per lipid. 4.2. Order Parameter. The deuterium order parameter provides a means to quantify ordering of the acyl tails, measuring the orientation of the deuterium−carbon (C−D) bonds with respect to the bilayer normal. Order parameters are of particular interest for validation of molecular mechanics force fields because they can be measured directly from deuterium NMR experiments,72,73 and unlike quantities such as area per lipid,74 the experimental data do not require models for interpretation.72 In molecular dynamics simulations, the order parameter is calculated via

Figure 5. Snapshot from NPT molecular dynamics simulations illustrating the complexation of DLPG with Na+. Upper image corresponds to the bilayer surface, while the lower image presents a side view. Atoms are colored as follows: Na+ (yellow), oxygen (red), phosphorus (gold), carbon (silver). Water molecules and free Na+ ions have been removed for clarity.

SCD =

9916

3 cos2 θ − 1 2

(5)



Article

Figure 6. Radial distribution functions for interaction of sodium ion with O2L (phosphate sp2 oxygen), OSLP (phosphate sp3 oxygen), OBL (glycerol sp2 oxygen), and PL (phosphorus) in DLPG: Roux Na+8 (black); Aqvist Na+20 (red).

Figure 8. Deuterium order parameters predicted by the TraPPE force field for DPPC, DMPC, DLPC, and POPC. Data are shown for the sn1 (open red squares) and sn-2 (open red circles) chains. Data from CHARMM C36 force field (green circles) and experimental data (solid black triangles).72,77,78

Table 2. Fraction of Bound Na+ and Lipids Predicted by Simulations of DLPG Performed with the TraPPE Force Field +

fraction of Na bound fraction of lipids bound no. of lipids/Na+

Na+ (Roux)

Na+ (Aqvist)

0.59 0.67 1.12

0.45 0.52 1.17

was observed between the predictions of the TraPPE force field and experimental data, with the TraPPE force field correctly reproducing all major features of the order parameter, such as the location of the kink in the sn-2 chain of POPC and C2 splitting in the sn-2 chain of all PC lipids. Compared to experimental data and the C36 force field, TraPPE predicts slightly less ordering, which may be due to limitations in the united-atom representation. Order parameter splitting at C2 position of sn-2 chain is a known feature for both saturated and unsaturated lipids. While not observed in older lipid force fields, such as the CHARMM C27 parameter set,37 nearly all modern all-atom force fields are capable of reproducing this phenomena.37,44,79 Reproduction of experimentally observed C2 splitting is the result of optimization of the torsional potentials governing dihedral rotations in the glycerol region. Similar splitting was observed for all PC-containing lipids. For DLPC, DMPC and DPPC, the lower splitting value was in excellent agreement with the experimental data, while the higher value is underpredicted by the TraPPE force field. For POPC, the magnitude of the splitting is reproduced, but ordering is underpredicted slightly. Experimentally, the SCD value for the sn-1 at C2 must be higher than for SCD for the sn-2 chain,73 and these results are also reproduced by the TraPPE force field. Order parameters for both sn-1 and sn-2 chains were calculated for DLPG, DMPG, DOPG, and DOPS and are shown in Figure 9. For DLPG and DMPG, simulations performed with the Aqvist Na+ parameters exhibit less ordering than similar calculations performed with the Roux Na+ parameters. This was expected, since the order parameter is a function of area per lipid, and the choice of Na+ ion parameters has a significant impact on the predicted area per lipid for DLPG and DMPG. The Na+ parameters have less of an effect on the predicted area per lipid for lipids that contain unsaturated chains, and the resulting order parameters for DOPG and DOPS display show similar insensitivity to the choice of Na+ parameters. 4.3. X-ray Scattering Form Factors and Electron Density. Like the deuterium order parameter, X-ray form factors F(q) may be determined directly from experimental data and provide an unambiguous measure of bilayer structure. In this work, structure factor for all lipids was determined using SIMtoEXP.80 Form factors predicted by simulation for DMPC,

Figure 7. Individual component electron density profiles for DLPG predicted by the TraPPE force field combined with the Roux Na+ (solid) and Aqvist Na+ (dashed) parameters: glycerol (orange), phosphate + choline (blue), methylene (red), terminal methyl (green), and water (black).

where θ is the angle formed by the C−D vector and the bilayer normal and the angular brackets denote an ensemble average.75 United-atom force fields, such as the ones developed in this work, lack hydrogen atoms to form the required Ci−Hi vectors. Therefore, hydrogen atoms were added to the lipid tails post simulation assuming ideal tetrahedral geometry. For saturated lipids, such as DPPC, equivalent results for the order parameter may be determined by using the vector joining Ci−1−Ci+1 carbons;75,76 the determination of order parameters for unsaturated lipid tails, however, requires the generation of pseudohydrogen atoms. Calculated values of order parameter are plotted as a function of the carbon number with their respective experimental values72,77,78 in Figure 8 for DPPC, DMPC, DLPC, and POPC. For DPPC and POPC, data from the CHARMM C36 force field are provided for comparison.37 Overall, close agreement 9917



Article

Figure 11. Form factors predicted by the TraPPE force field combined with Aqvist Na+ (black solid line) and Roux Na+ (green solid line) for DLPG, DMPG, DOPG, and DOPS; Experimental data for unilamellar and oriented vesicles are represented by red circles.

Figure 9. Order parameters as predicted by the TraPPE force field for DLPG, DMPG, DOPG, and DOPS for simulations using the Roux Na+ (black) and Aqvist Na+ (red) parameters. Sn-1chains are shown as open squares and sn-2 chains as open circles.

DPPC, POPC, and DLPC are shown in Figure 10. In each case, the TraPPE force field predicts peak and minima locations that

Figure 10. Form factors predicted by the TraPPE force field for DLPC, DMPC, DPPC, and POPC (black solid line) and compared with experimental data for unilamellar and oriented vesicles (red circles).

Figure 12. (Top) Total electron density profile for DPPC predicted by the TraPPE force field. Experimental data (black), simulation (red). (Bottom): Individual component density profiles for DPPC predicted by the TraPPE force field. Simulation (dashed) and experiment (solid) are shown for each component: glycerol (orange), phosphate + choline (blue), methylene (red), terminal methyl (green), and water (black).

are in close agreement with experimental data,81,82 while some discrepancies between simulation and experiment were observed for data at low q factors (q < 0.2 Å−1). This may be due to limitations of the united-atom representation, which displays slightly less ordering of the alkyl tails compared to allatom simulations and experimental data. X-ray form factors were determined for DLPG, DMPG, DOPG, and DOPS using the Roux and Aqvist Na+ parameters (Figure 11). In comparison to experimental data, both sets of ion parameters predict the locations of the peaks and minima in agreement with experimental data, although there are slight deviations. For all of the PG-containing lipids, the simulation data are shifted slightly to larger q values for the first peak. This shift is most pronounced for DMPG, and less so for DLPG and DOPG. A similar phenomena was observed in simulations of the Slipids force field,68 and it was suggested that this may be due to an underprediction of the area per lipid.83 In this case, however, the TraPPE force field predicts the area per lipid of DMPG to within 0.3 Å2 of experimental data, suggesting other factors may be responsible for the observed deviations between simulation and experiment. In Figure 12 the total electron density and the component electron densities are presented for DPPC and compared to

experimental data. Electron density profiles for other lipids are presented in the Supporting Information (Figures S12 and S13). Electron densities were extracted from the simulation data for each lipid using SIMtoEXP.80 It should be noted that transformation of experimental F(q) data into real space ρ(z) (electron density) requires the use of structural models.74,84,85 Close agreement between simulation and experiment was observed for the total electron density in the headgroup region, while subtle differences were observed near the center of the bilayer. Examination of the individual components of the electron density shows excellent agreement between simulation and experiment for the phosphate + choline group (PCN), glycerol region (CG), and water. The electron density for the terminal CH3 group is substantially narrower than experiment. The largest difference between simulation and experiment was observed for the CH2 group, where the simulated electron density is overpredicted and the trough at the bilayer center is not as deep. This is another manifestation of the united-atom approximation. 9918



Article

For all the lipids, the bilayer thickness DHH was determined by using the peak to peak distance given by the total electron density. These data are also listed in Table 1. For the PC lipids, an average absolute deviation (AAD) of 0.42 Å was observed, which is an improvement over the GROMOS 53A6 parameter set (AAD = 2.12 Å), and compares favorably with the all-atom force fields Slipids (AAD = 0.48 Å) and GAFFLipid (AAD = 0.88 Å). For PG lipids, TraPPE + Aqvist Na+ predicts an AAD of 0.23 Å vs an AAD of 1.3 Å reported for the Slipids force field.68

optimization of partial charges and bonded parameters, such as barriers to dihedral rotation.

■

ASSOCIATED CONTENT

* Supporting Information S

Additional data for rotational barriers, evolution of volume per lipid with simulation time, and electron densities may be found in the Supporting Information. The force field parameters for the various models developed in this study are supplied in a file format compatible with NAMD. This material is available free of charge via the Internet at http://pubs.acs.org..

■

5. CONCLUSIONS In this work, the Transferable Potentials for Phase Equilibria (TraPPE) force field was extended to phospholipids. Models were presented for a variety of lipid head groups (PC, PE, PG, and PS) and saturated and unsaturated tails containing 12−18 carbon atoms. The predictions of the TraPPE force field for the area per lipid were within 1−2% of experimental values except for POPE (5.6% deviation), and POPC and POPS (3% deviation). Similar predictive capabilities were observed for quantities such as bilayer thickness and volume per lipid. Key structural properties of the bilayer such as order parameter splitting in the sn-2 chain and X-ray form factors were also reproduced closely. Additional calculations were performed to elucidate the role of water and ion models on the performance of the phospholipid force field. Simulations with the TIP3P force field showed deeper penetration of water in to the lipid bilayer, leading to increased values of the area per lipid compared to the SPC/E water model. The choice of Na+ parameters was also found to influence the hydration of the lipid head groups. Simulations performed with the Roux parametrization displayed greater binding to oxygen atoms in the headgroup compared to simulations performed with parameters from Aqvist. Increased binding of Na+ to lipid head groups produced a slight dehydration effect, lowering the observed area per lipid and increasing ordering of the alkyl tails. Overall, these results highlight the broad transferability of the Transferable Potentials for Phase Equilibria, which was originally developed for the prediction of vapor−liquid coexistence, to the simulation of complex biological systems. The accurate prediction of vapor−liquid coexistence requires optimization of nonbonded parameters to produce accurate results over large regions of pressure−volume−temperature space. In many cases, parameters are optimized to reproduce experimental data over a 200 K temperature range and a 50 bar range in pressure. Changes of 0.01 kcal/mol in the LennardJones parameters produce measurable changes in the predicted phase diagram. Therefore, it is expected that nonbonded parameters optimized to reproduce vapor−liquid coexistence will exhibit significantly better transferability to other conditions than parameters optimized to reproduce physical properties at a single temperature and pressure (usually 298 K and 1.013 bar). It is interesting to note that a number of recent revisions to lipid force fields, such as AMBER34 and GROMOS 43A1-S3,25 and Slipids,44 have focused on a reparametrization of the lipid tails to reproduce liquid densities and heats of vaporization for longer alkanes, such as pentadecane. While nonbonded parameters of the TraPPE force field have already undergone rigorous optimization, the further application of TraPPE to biological molecules, such as nucleic acids and proteins, will require significant effort focused on the

AUTHOR INFORMATION

Corresponding Author

*E-mail: jpotoff@wayne.edu. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank John Nagle and Norbert Kučerka for providing x-ray form factor data, and providing technical assistance on the use of SIMtoEXP for the extraction of X-ray structure factors from the simulation data. Funding from the National Science Foundation NSF CBET-1066661, a Thomas C. Rumble Fellowship (N.B.), and the Wayne State University Research Enhancement Program are gratefully acknowledged. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant no. OCI-1053575. Additional computational resources were provided by the Grid computational environment at Wayne State University.

■

REFERENCES

(1) Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J. J. Am. Chem. Soc. 1996, 118 (45), 11225−11236. (2) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117 (19), 5179−5197. (3) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102 (18), 3586−3616. (4) Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. J. Comput. Chem. 2004, 25 (13), 1656−1676. (5) Halgren, T. A. J. Comput. Chem. 1995, 17, 520−552. (6) Yin, D. X.; Mackerell, A. D. J. Comput. Chem. 1998, 19 (3), 334− 348. (7) Chen, I. J.; Yin, D. X.; MacKerell, A. D. J. Comput. Chem. 2002, 23 (2), 199−213. (8) Beglov, D.; Roux, B. J. Chem. Phys. 1994, 100 (12), 9050−9063. (9) Garrido, N. M.; Jorge, M.; Queimada, A. J.; Gomes, J. R. B.; Economou, I. G.; Macedo, E. A. Phys. Chem. Chem. Phys. 2011, 13 (38), 17384−17394. (10) Xu, Z. T.; Luo, H. H.; Tieleman, D. P. J. Comput. Chem. 2007, 28 (3), 689−697. (11) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102 (14), 2569−2577. (12) Nath, S. K.; Escobedo, F. A.; de Pablo, J. J. J. Chem. Phys. 1998, 108 (23), 9905−9911. (13) Errington, J. R.; Panagiotopoulos, A. Z. J. Phys. Chem. B 1999, 103 (30), 6314−6322. 9919



Article

(49) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2004, 108 (45), 17596−17605. (50) Wick, C. D.; Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 2000, 104 (33), 8008−8016. (51) Wick, C. D.; Stubbs, J. M.; Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2005, 109 (40), 18974−18982. (52) Sokkalingam, N.; Kamath, G.; Coscione, M.; Potoff, J. J. J. Phys. Chem. B 2009, 113 (30), 10292−10297. (53) Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2007, 111 (36), 10790− 10799. (54) Kamath, G.; Cao, F.; Potoff, J. J. J. Phys. Chem. B 2004, 108 (37), 14130−14136. (55) Hogberg, C. J.; Nikitin, A. M.; Lyubartsev, A. P. J. Comput. Chem. 2008, 29 (14), 2359−2369. (56) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11 (3), 361−373. (57) 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.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision A.02; Gaussian Inc: Wallingford, CT, 2009. (58) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91 (24), 6269−6271. (59) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. J. Comput. Chem. 2005, 26 (16), 1781−1802. (60) Bhatnagar, N.; Kamath, G.; Chelst, I.; Potoff, J. J. J. Chem. Phys. 2012, 137 (1), 014502. (61) Martyna, G. J.; Tobias, D. J.; Klein, M. L. J. Chem. Phys. 1994, 101 (5), 4177−4189. (62) Feller, S. E.; Zhang, Y. H.; Pastor, R. W.; Brooks, B. R. J. Chem. Phys. 1995, 103 (11), 4613−4621. (63) Taylor, J.; Whiteford, N. E.; Bradley, G.; Watson, G. W. Biochim. Biophys. Acta, Biomembr. 2009, 1788 (3), 638−649. (64) Rappolt, M.; Hickel, A.; Bringezu, F.; Lohner, K. Biophys. J. 2003, 84 (5), 3111−3122. (65) Cordomi, A.; Edholm, O.; Perez, J. J. J. Chem. Theory Comput. 2009, 5 (8), 2125−2134. (66) Bockmann, R. A.; Hac, A.; Heimburg, T.; Grubmuller, H. Biophys. J. 2003, 85 (3), 1647−1655. (67) Gurtovenko, A. A.; Vattulainen, I. J. Phys. Chem. B 2008, 112 (7), 1953−1962. (68) Jambeck, J. P. M.; Lyubartsev, A. P. J. Chem. Theory Comput. 2013, 9 (1), 774−784. (69) Poger, D.; Mark, A. E. J. Chem. Theory Comput. 2010, 6 (1), 325−336. (70) Feller, S. E.; Pastor, R. W. Biophys. J. 1996, 71 (3), 1350−1355. (71) Waheed, Q.; Edholm, O. Biophys. J. 2009, 97 (10), 2754−2760. (72) Petrache, H. I.; Dodd, S. W.; Brown, M. F. Biophys. J. 2000, 79 (6), 3172−3192. (73) Seelig, A.; Seelig, J. Biochemistry 1974, 13 (23), 4839−4845. (74) Kucerka, N.; Nagle, J. F.; Sachs, J. N.; Feller, S. E.; Pencer, J.; Jackson, A.; Katsaras, J. Biophys. J. 2008, 95 (5), 2356−2367. (75) Vermeer, L. S.; de Groot, B. L.; Reat, V.; Milon, A.; Czaplicki, J. Eur. Biophys. J. 2007, 36 (8), 919−931. (76) Ulmschneider, J. P.; Ulmschneider, M. B. J. Chem. Theory 2009, 5 (7), 1803−1813.

(14) Ungerer, P.; Beauvais, C.; Delhommelle, J.; Boutin, A.; Rousseau, B.; Fuchs, A. H. J. Chem. Phys. 2000, 112 (12), 5499−5510. (15) Potoff, J. J.; Bernard-Brunel, D. A. J. Phys. Chem. B 2009, 113 (44), 14725−14731. (16) Kamath, G.; Lubna, N.; Potoff, J. J., J. Chem. Phys. 2005, 123 (12). (17) Ketko, M. H.; Kamath, G.; Potoff, J. J. J. Phys. Chem. B 2011, 115 (17), 4949−4954. (18) Potoff, J. J.; Siepmann, J. I. AIChE J. 2001, 47 (7), 1676−1682. (19) Kamath, G.; Georgiev, G.; Potoff, J. J. J. Phys. Chem. B 2005, 109 (41), 19463−19473. (20) Aqvist, J. J. Phys. Chem. 1990, 94 (21), 8021−8024. (21) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25 (9), 1157−1174. (22) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; MacKerell, A. D. J. Comput. Chem. 2010, 31 (4), 671−690. (23) Berger, O.; Edholm, O.; Jahnig, F. Biophys. J. 1997, 72 (5), 2002−2013. (24) Chiu, S. W.; Clark, M.; Balaji, V.; Subramaniam, S.; Scott, H. L.; Jakobsson, E. Biophys. J. 1995, 69 (4), 1230−1245. (25) Chiu, S. W.; Pandit, S. A.; Scott, H. L.; Jakobsson, E. J. Phys. Chem. B 2009, 113 (9), 2748−2763. (26) Poger, D.; Van Gunsteren, W. F.; Mark, A. E. J. Comput. Chem. 2010, 31 (6), 1117−1125. (27) Dickey, A.; Faller, R. Biophys. J. 2008, 95 (6), 2636−2646. (28) Kukol, A. J. Chem. Theory Comput. 2009, 5 (3), 615−626. (29) Zhao, W.; Rog, T.; Gurtovenko, A. A.; Vattulainen, I.; Karttunen, M. Biophys. J. 2007, 92 (4), 1114−1124. (30) Piggot, T. J.; Holdbrook, D. A.; Khalid, S. J. Phys. Chem. B 2011, 115 (45), 13381−13388. (31) Cieplak, P.; Cornell, W. D.; Bayly, C.; Kollman, P. A. J. Comput. Chem. 1995, 16 (11), 1357−1377. (32) Jojart, B.; Martinek, T. A. J. Comput. Chem. 2007, 28 (12), 2051−2058. (33) Skjevik, A. A.; Madej, B. D.; Walker, R. C.; Teigen, K. J. Phys. Chem. B 2012, 116 (36), 11124−11136. (34) Dickson, C. J.; Rosso, L.; Betz, R. M.; Walker, R. C.; Gould, I. R. Soft Matter 2012, 8 (37), 9617−9627. (35) Feller, S. E.; Venable, R. M.; Pastor, R. W. Langmuir 1997, 13 (24), 6555−6561. (36) Husslein, T.; Newns, D. M.; Pattnaik, P. C.; Zhong, Q. F.; Moore, P. B.; Klein, M. L. J. Chem. Phys. 1998, 109 (7), 2826−2832. (37) 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. J. Phys. Chem. B 2010, 114 (23), 7830−7843. (38) Feller, S. E.; Pastor, R. W. J. Chem. Phys. 1999, 111 (3), 1281− 1287. (39) Feller, S. E.; MacKerell, A. D. J. Phys. Chem. B 2000, 104 (31), 7510−7515. (40) Jensen, M. O.; Mouritsen, O. G.; Peters, G. H. Biophys. J. 2004, 86 (6), 3556−3575. (41) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97 (40), 10269−10280. (42) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollman, P. A. J. Am. Chem. Soc. 1993, 115 (21), 9620−9631. (43) Sonne, J.; Jensen, M. O.; Hansen, F. Y.; Hemmingsen, L.; Peters, G. H. Biophys. J. 2007, 92 (12), 4157−4167. (44) Jambeck, J. P. M.; Lyubartsev, A. P. J. Phys. Chem. B 2012, 116 (10), 3164−3179. (45) Jambeck, J. P. M.; Lyubartsev, A. P. J. Chem. Theory Comput. 2012, 8 (8), 2938−2948. (46) Chen, B.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2001, 105 (15), 3093−3104. (47) Lubna, N.; Kamath, G.; Potoff, J. J.; Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2005, 109 (50), 24100−24107. (48) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1999, 103 (21), 4508−4517. 9920



Article

(77) Perly, B.; Smith, I. C. P.; Jarrell, H. C. Biochemistry 1985, 24 (4), 1055−1063. (78) Seelig, J.; Waespesarcevic, N. Biochemistry 1978, 17 (16), 3310− 3315. (79) Siu, S. W. I.; Vacha, R.; Jungwirth, P.; Bockmann, R. A. J. Chem. Phys. 2008, 128 (12). (80) Kucerka, N.; Katsaras, J.; Nagle, J. F. J. Membr. Biol. 2010, 235 (1), 43−50. (81) Petrache, H. I.; Tristram-Nagle, S.; Gawrisch, K.; Harries, D.; Parsegian, V. A.; Nagle, J. F. Biophys. J. 2004, 86 (3), 1574−1586. (82) Kucerka, N.; Liu, Y. F.; Chu, N. J.; Petrache, H. I.; TristramNagle, S. T.; Nagle, J. F. Biophys. J. 2005, 88 (4), 2626−2637. (83) Kucerka, N.; Holland, B. W.; Gray, C. G.; Tomberli, B.; Katsaras, J. J. Phys. Chem. B 2012, 116 (1), 232−239. (84) Klauda, J. B.; Kucerka, N.; Brooks, B. R.; Pastor, R. W.; Nagle, J. F. Biophys. J. 2006, 90 (8), 2796−2807. (85) Sachs, J. N.; Petrache, H. I.; Woolf, T. B. Chem. Phys. Lipids 2003, 126 (2), 211−223. (86) Nagle, J. F.; Tristram-Nagle, S. Biochim. Biophys. Acta, Rev. Biomembr. 2000, 1469 (3), 159−195. (87) Poger, D.; Mark, A. E. J. Chem. Theory Comput. 2012, 8 (11), 4807−4817. (88) Rawicz, W.; Olbrich, K. C.; McIntosh, T.; Needham, D.; Evans, E. Biophys. J. 2000, 79 (1), 328−339. (89) Pan, J.; Tristram-Nagle, S.; Kucerka, N.; Nagle, J. F. Biophys. J. 2008, 94 (1), 117−124. (90) Kucerka, N.; Tristram-Nagle, S.; Nagle, J. F. J. Membr. Biol. 2005, 208 (3), 193−202. (91) Binder, H.; Gawrisch, K. J. Phys. Chem. B 2001, 105 (49), 12378−12390. (92) Rand, R. P.; Parsegian, V. A. Biochim. Biophys. Acta 1989, 988 (3), 351−376. (93) Pan, J. J.; Heberle, F. A.; Tristram-Nagle, S.; Szymanski, M.; Koepfinger, M.; Katsaras, J.; Kucerka, N. BBA-Biomembranes 2012, 1818 (9), 2135−2148. (94) Henin, J.; Shinoda, W.; Klein, M. L. J. Phys. Chem. B 2009, 113 (19), 6958−6963.

9921


Biomolecular Simulations with the Transferable Potentials for Phase

Recommend Documents