CHARMM36 United Atom Chain Model for Lipids and Surfactants

Article pubs.acs.org/JPCB

CHARMM36 United Atom Chain Model for Lipids and Surfactants Sarah Lee,†,∥ Alan Tran,†,∥ Matthew Allsopp,†,∥ Joseph B. Lim,† Jérôme Hénin,*,‡ and Jeffery B. Klauda*,† †

Department of Chemical and Biomolecular Engineering, University of Maryland, College Park, Maryland 20742, United States Laboratoire de Biochimie Théorique, CNRS, Institut de Biologie Physico-Chimique, 13 rue Pierre et Marie Curie, 75005 Paris, France

‡

S Supporting Information *

ABSTRACT: Molecular simulations of lipids and surfactants require accurate parameters to reproduce and predict experimental properties. Previously, a united atom (UA) chain model was developed for the CHARMM27/27r lipids (Hénin, J., et al. J. Phys. Chem. B. 2008, 112, 7008−7015) but suffers from the flaw that bilayer simulations using the model require an imposed surface area ensemble, which limits its use to pure bilayer systems. A UA-chain model has been developed based on the CHARMM36 (C36) all-atom lipid parameters, termed C36-UA, and agreed well with bulk, lipid membrane, and micelle formation of a surfactant. Molecular dynamics (MD) simulations of alkanes (heptane and pentadecane) were used to test the validity of C36-UA on density, heat of vaporization, and liquid self-diffusion constants. Then, simulations using C36-UA resulted in accurate properties (surface area per lipid, X-ray and neutron form factors, and chain order parameters) of various saturated- and unsaturated-chain bilayers. When mixed with the all-atom cholesterol model and tested with a series of 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC)/cholesterol mixtures, the C36-UA model performed well. Simulations of self-assembly of a surfactant (dodecylphosphocholine, DPC) using C36-UA suggest an aggregation number of 53 ± 11 DPC molecules at 0.45 M of DPC, which agrees well with experimental estimates. Therefore, the C36-UA force field offers a useful alternative to the all-atom C36 lipid force field by requiring less computational cost while still maintaining the same level of accuracy, which may prove useful for large systems with proteins. CHARMM36 (C36)4 force field. In this force field, only polar hydrogens on the hydroxyl or NH are explicitly included. However, it has been shown that molecular simulations using the Berger-based force field7 yield bilayers too stable for ceramide lipids.13 Hénin et al.6 developed a UA acyl-chain model for the C27/ C27r14−17 lipids (C27-UA). This model maintained the important AA CHARMM character of the lipid headgroup while reducing the AA description of the chains. Because this was based on the C27r lipids, it required MD simulations in the constant area ensemble as this older version of the CHARMM lipids did not accurately represent lateral density and phase behavior of certain lipids in the constant number of molecules, temperature, and pressure (NPT) ensemble.4,18,19 This limits the ability to use C27-UA for lipid mixtures and lipids with proteins for which the average surface area (SA) per lipid is unknown. Presented in this paper is a redeveloped UA-chain parameter set that is compatible with C36 and allows NPT MD simulations, referred to here as C36-UA. C36-UA is an alternative to the AA C36 force field and has been parametrized such that it is compatible with the AA CHARMM family of

1. INTRODUCTION Lipids are an important building block in biology that typically form lamellar structures surrounding cells and organelles.1 Depending on its composition, a lipid mixture can form liquiddisordered (Ld) or liquid-ordered (Lo) bilayer phases; certain mixtures with sterols can result in Ld/Lo coexistence with domains of Lo.1,2 Molecular dynamics (MD) simulations have been used to study properties of lipid membranes, as they describe and predict structural detail inaccessible to experimental measurements. Essential to MD is an accurate description of how molecules interact via a mathematical representation of forces, known as a force field.3 These molecular-level force fields can be classified into three groups, i.e., all-atom (AA, including hydrogen),4,5 united atom (UA, hydrogens grouped with heavy atoms),6−9 and coarse grained (CG, grouping of heavy atoms).10−12 AA force fields require the most computational cost but provide the most detail compared to UA and CG force fields. AA force fields are limited to at best microsecond time scales on traditional supercomputers, whereas CG force fields can reach the ms time scale. GROMOS-based force fields have been the most used and accurate UA model for lipids.7−9 This UA grouping results in a reduction of interaction sites and computational cost, allowing MD simulations using the UA model to run for a longer time scale compared to what is possible using the all-atom © 2013 American Chemical Society

Received: October 18, 2013 Revised: December 13, 2013 Published: December 16, 2013 547

dx.doi.org/10.1021/jp410344g | J. Phys. Chem. B 2014, 118, 547−556

The Journal of Physical Chemistry B

Article

Table 1. Details of the Bilayer Simulations no. of atoms no. of lipids no. of waters temp. (K) prod. time (ns) a

DMPC

DPPC

POPC

DOPC

DMPC/chola

10 440 72 1848 303.0 50

11 751 72 2261 323.0 50

12 054 72 2242 303.0 50

12 699 72 2409 303.0 50

20 874, 15 854, 15 974 132, 100, 100 3958, 2998, 2998 308.15, 308.15, 298.15 50, 50, 100

For mixtures with cholesterol (chol), three concentrations were simulated and are reported in the order 3, 10, and 30% cholesterol.

Table 2. Details of the DPC Micelle Simulationsa no. of atoms no. of DPCs no. of waters [DPC]b prod. time (ns) formation time of micelle(s) (ns)

M1

M2

M3

M4

M5

M6

22 053 54 6 667 0.45 150 56,c 88

44 103 108 13 333 0.45 150 41

66 156 162 20 000 0.45 150 81

34 093 20 11 111 0.1 150 39

68 186 40 22 222 0.1 150 88d

27 684 36 8 772 0.23 150 52

a

M1−M6 refer to the micelle simulation. bThe concentration of DPC ([DPC]) is in moles per liter. cM1 micelle fusion. dM5 never equilibrates, but at this time, the two micelles have formed.

force fields. First, the methods used to develop the UA force field parameters will be described in Section 2. The methods used to run and analyze the MD simulations of bulk alkanes, lipid bilayers, and micelles are also described in this section. Results will be presented in Section 3, starting with optimization of the UA parameters in 3.1. Then, results from single-lipid bilayers and mixed bilayers with AA cholesterol will be presented afterward. The results will conclude with MD simulations of the self-assembly of a detergent into micelles under various conditions. These results will be discussed in Section 4 focusing on the accuracy of C36-UA, insight into comparisons with X-ray/neutron scattering experiments, and micelle formation.

conditions as the lipid simulations described below, at 1 bar and 312 K for pentadecane and 310 K for the other compounds. Torsion angle histograms were collected. All MD simulations were performed using NAMD21 with 2 fs timesteps. A force-based switching function was used to switch the van der Waals forces to zero over a range of 8−12 Å. Longrange electrostatic interactions were included with the particlemesh Ewald method22 with an interpolation order of 4 and a direct space tolerance of 10−6. Langevin dynamics was used to maintain constant temperatures for each system, while the Nosé−Hoover Langevin-piston algorithm23,24 was used to maintain constant pressure at 1 bar. For the isotropic simulations (alkanes and micelles), the unit cell was cubic with X = Y = Z. For simulations with lipid bilayers, the semiisotropic cell was used with X constrained to be equal to Y but allowed to vary independently with respect to Z. The alkane simulations were run with 128 and 64 heptane and pentadecane molecules, respectively. Details regarding the setup of the lipid bilayers and micelle simulations are given in Tables 1 and 2, respectively. MD simulations of mixtures with AA cholesterol used the updated cholesterol force field parameters from Lim et al.25 The C36 lipid force field4 was used for the non-UA portion of the lipids. The modified TIP3P water model26,27 was used for water to be consistent with the CHARMM family of force fields. Previously equilibrated coordinates for alkanes6 were used as starting coordinates in our simulations. The energy was minimized with the conjugate gradient method in NAMD for 200 steps. Thermal equilibration was performed for 1 ns and productions runs were each run for 10 ns. For the bilayers, previously equilibrated AA coordinates4,6,25 were used as starting structures for the C36-UA simulations. For the DMPC/cholesterol systems, the AA bilayers were initially built with CHARMM-GUI28,29 and equilibrated before converting to the UA topology. The lipids without cholesterol were energy minimized for 200 steps and those with cholesterol for 500 steps. Thermal equilibration was performed for 1 ns, and the lengths of the production runs are reported in Table 1. In addition to bilayer simulations, micelle simulations of ndodecylphosphocholine (DPC) were run. A random mixture of DPC in water was set up with PACKMOL30 and used to

2. METHODS A set of four biobased UA force fields (OPLS,20 Berger,7 and Chiu et al. published in 19998 and 20039) were used to evaluate their ability to match the experimental data (densities, heats of vaporization, and diffusion) of alkanes, i.e., heptane and pentadecane. Although these force fields were tested with the C27-UA6 development, C36 requires the use of a force-based switching function that alters the good agreement of UA parameters for alkanes developed by Chiu et al.9 To improve the description in the critical ester region, explicit hydrogen atoms were reintroduced in the first methylene group of each acyl chain. This restores the carbon/hydrogen charge separation on this group and places the link between the AA and UA regions at a bond between methylene groups. Bonded interactions between the AA and UA regions at the seam are described by standard C36 parameters for bonds and angles and specific torsional potentials that are key to reproducing the AA conformational preferences. As with C27-UA, the torsional terms were fit using potentials of mean force calculations.6 Torsional parameters were optimized iteratively to reproduce torsional free energy profiles for several model compounds: pentadecane for saturated hydrocarbon chains, cis-5-decene for monounsaturated chains, and methyl hexanoate for the ester region. The functional form of C27-UA torsions, order-6 cosine series, was retained. Torsional distributions were computed from bulk liquid simulations of the pure model compounds, in the same 548



Article

Table 3. Bulk Properties of Alkane Simulations at 312.15 K and 1 bara OPLS20

Berger7

C-19998

C-20039

expt

n-heptane ρ ΔHvap DPBC Ds

0.64 7.21 2.15 2.76

± ± ± ±

0.00 0.11 0.01 0.01

0.52 4.60 3.21 3.78

± ± ± ±

0.00 0.10 0.02 0.02

ρ ΔHvap DPBC Ds

0.77 16.5 0.59 0.70

± ± ± ±

0.00 0.1 0.01 0.01

0.69 11.5 0.90 1.01

± ± ± ±

0.00 0.2 0.01 0.01

0.59 4.91 3.11 3.70

± ± ± ±

0.00 0.08 0.01 0.01

0.66 6.80 2.28 2.89

± ± ± ±

0.00 0.11 0.01 0.01

0.6736 8.5037

0.72 11.8 0.94 1.04

± ± ± ±

0.00 0.1 0.00 0.00

0.73 14.8 0.68 0.79

± ± ± ±

0.00 0.2 0.00 0.00

0.7636 17.837

3.6838

n-pentadecane

−5

0.66b

Density is in grams per cubic centimeter. ΔH is in kilocalories per mole. Diffusion is in 10 cm /s. ± standard errors are reported (values of zero are smaller than reported decimal place). bDiffusion estimated at 312.15 K from given experimental temperature range of 288−303 K.35 a

vap

methylene hydrogens on the second carbon. Complete topology and parameter files are available online from ref 63 and are also included in the Supporting Materials. 3.2. Single-Lipid Bilayers. The SA/lipid is a common metric to determine equilibration in the lateral density of lipid bilayers. The four phospholipid bilayers all resulted in equilibrated values for SA/lipid within 10 ns (Figure S1 of the Supporting Information). The averages for C36-UA agreed with the AA C36 and X-ray and neutron scattering (Table 4).

simulate micelle formation. Initially, 200 steps of energy minimization were used, followed by 1 ns of thermal equilibration at 300 K. MD simulations were run in a cubic simulation cell for a production length of 150 ns each. More details regarding the MD micelle simulations are reported in Table 2. For the alkane simulations, bulk properties of density, heat of vaporization, and self-diffusion were calculated. The selfdiffusion (Ds) was estimated with the following equation to include system-size effects using a hydrodynamic model31 Ds = DPBC +

kBTξ 6πηL

2

Table 4. Average Surface Area per Lipid (Å2) for Bilayers with ± Standard Errors

(1)

C36-UA

where DPBC is the uncorrected diffusion constant, L the box length, and ξ = 2.837297. Because the bilayer system sizes were small, undulations can be neglected32 and the surface area per lipid (SA/lipid) was just taken as the total X·Y divided by the number of lipids per leaflet. All averages for bilayer properties were taken from the last 25 ns, unless otherwise specified. The deuterium order parameters were estimated by reconstructing the deuterium atoms assuming an ideal geometry. The SIMtoEXP software33 was used to obtain the X-ray and neutron scattering form factors (F(q)) and electron density profiles (EDPs). The Visual Molecular Dynamics (VMD)34 program was used to create snapshots.

DMPC DPPC POPC DOPC

61.6 61.7 64.4 67.5

± ± ± ±

0.3 0.3 0.3 0.2

C364 60.8 62.9 64.7 69.0

± ± ± ±

0.2 0.3 0.2 0.3

expt 60.6 63.0 68.3 67.4

± ± ± ±

0.540 1.039 1.541 1.039

The differences between C36 and C36-UA were negligible for 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) and 1palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC). MD simulations of the shorter-chained DMPC lipid (14 carbons versus 16 in DPPC) resulted in a higher SA/lipid by 0.8 Å2 compared to that of C36. For 1,2-dioleoyl-sn-glycero-3phosphocholine (DOPC), the SA/lipid decreased by 1.5 Å2 compared to that of C36 but was in excellent agreement with a model fit to scattering experiments.39 Because the experimentally based SA/lipid results from a structural model fit to the raw experimental data, deficiencies in the structural model may influence the accuracy. Therefore, a more direct comparison is with the X-ray and neutron scattering F(q).18,42,43 Comparisons are shown in Figure 1 for DOPC and for the other lipid bilayers in Figures S2−S4 of the Supporting Information. Overall the agreement between experiment and C36-UA was satisfactory. The largest errors corresponded with DMPC and were consistent with the SA/lipid being too high for DMPC, i.e., the crossing point for the X-ray F(q) was larger for C36-UA compared to that of experiment (Figure S2 of the Supporting Information). Agreement between C36-UA and both X-ray and neutron scattering suggests a well-parametrized UA force field in terms of bilayer profile density. This model is somewhat worse than C36 for DMPC but an improvement for DOPC. The EDPs are in Figure 2 for DOPC and Figure S5 of the Supporting Information for the other lipids. The chain density for C36-UA was consistently higher than that of C36, which is

3. RESULTS 3.1. Optimizing C36-UA Parameters. Results from MD simulations using the four UA force fields are shown in Table 3. On the basis of the experimental results, OPLS was chosen as the best parameter to use for lipid bilayer simulations. OPLS generated accurate density values for n-heptane and npentadecane when compared to the literature values, with errors of 3.0 and 1.3%, respectively. The heat of vaporization (ΔHvap) was consistently lower than the experimental values for all UA force fields. However, OPLS resulted in the closest estimate for ΔH vap with 15.1% (heptane) and 7.3% (pentadecane) error from experiment. Because lipids tend to have acyl-chain lengths of 12−18 carbons, the reasonable agreement with ΔHvap was deemed acceptable. Although the Ds was slower than experiment for n-heptane, the agreement for npentadecane was excellent. Therefore, C36-UA uses the OPLS values for the 6−12 Lennard-Jones potential. Refinement of the torsional parameters from C27-UA entailed only small changes in the alkane parameters (at most 0.05 kcal/mol). Parameters involving the ester side of the chain had to be rederived entirely because of the addition of two 549



Article

due to the UA potential used. This comes from the OPLS CH2 group (Figure S5 of the Supporting Information) that has a density higher than that of C36 and resulted in a lower component volume (25 versus 28 Å3/group). The total volume/lipid was also smaller than the experimental volumes,44 accounting for the discrepancy with the value of the X-ray F(0) in Figure 1 using the formula in ref 45. Aside from this difference, C36 and C36-UA agreed well. The DMPC head-toheadgroup spacing (DHH) was smaller for C36-UA, consistent with the larger SA/lipid. The final verification of the C36-UA force field was comparison of chain order parameters (Figure 3). The

Figure 1. DOPC form factors from X-ray (top) and neutron (bottom) scattering.42,43 ULV = unilamellar vesicles. F(q = 0) = 0 for X-ray.

Figure 3. SCD’s for DPPC46,47 and POPC50,51 from experiment are shown in red symbols (sn-1, triangles; sn-2, squares). Blue is for MD simulations with sn-1 chain and green for sn-2. C364 is shown in dashed lines, while C36-UA is in solid. Error bars are similar to symbol sizes and not shown.

headgroup is still an AA potential in C36-UA, so only the chain SCDs were compared. C36 was optimized to obtain an accurate representation of DPPC SCDs and C36-UA agreed similarly with experiment.46,47 Overall, the C36-UA SCDs were slightly lower than those of C36 (Figure 3). C36-UA results indicated that the sn-2 chain had a lower SCD for C2 compared to that of sn-1, which is in agreement with experiment.47 Comparing experiments for sn-1 and sn-2 deuterated separately,48,49 sn-2 had a higher order for higher carbons down the chain and MD simulations using C36-UA agreed with experiment. POPC SCDs for C36-UA agreed well with experiment50,51 and were of similar accuracy compared to those for C36 (Figure 3, bottom panel). The sn-2 chain was better represented by C36-UA. For the sn-1 chain, C36-UA matched experiment below C9 but underestimated the order for carbons closer to the headgroup. Overall, the dihedral fits for the UA potential appeared to accurately represent the chain order for saturated and unsaturated PC lipids. 3.3. Phospholipids with Cholesterol. With a wellparametrized chain UA model, our next step to validate this

Figure 2. Symmetrized EDP for the SDP model39 based on fits to F(q) is shown in red for DOPC. EDPs from MD simulations are shown in black for C36-UA and C36.4

550



Article

model was to determine how well it works with the AA sterol force field.25 The time series for the SA/lipid for the three cholesterol concentrations with DMPC are shown in Figure S6 of the Supporting Information. The two lowest concentrations equilibrated quickly, and averages were calculated from the last 25 ns. For DMPC/30% cholesterol, the SA/lipid slowly decreased until about 75 ns, and data from the last 25 ns of the 100 ns run was used for analysis. The average SA/lipid was 60.5 ± 0.2, 55.5 ± 0.2, and 41.4 ± 0.1 Å2 for 3, 10, and 30% cholesterol, respectively. The SA/lipid values were similar (slightly lower) to those based on our previous AA C36 force field with its updated cholesterol parameters, i.e., 56.1 ± 0.1 and 43.0 ± 0.1 Å2 for 10 and 30% cholesterol, respectively. The X-ray form factors from C36-UA compared favorably with experiments52 for 3 and 10% cholesterol (Figure S7 of the Supporting Information). The peak positions were shifted to the right for these bilayers, suggesting that the SA/lipid may be too high in the C36-UA force field (similar to DMPC-only bilayers as shown in Figure S2 of the Supporting Information). The real-space electron densities for the 10 and 30% cholesterol bilayers are shown in Figure 4. As with pure DMPC (Figure S5

Figure 5. SCD’s for DMPC/cholesterol bilayers for C36-UA, C36,25 and experiment for the sn-2 chain.48 Error bars are similar to symbol sizes and not shown.

Similar to our previous simulations, the lipid order was lower than experiment for DMPC/30% cholesterol. The lower SA/ lipid for C36-UA did not appear to influence the SCDs at this sterol concentration, which may be the result of the DMPC lipid taking up less volume. This lower volume would result in a lower SA/lipid and is consistent with the higher density in the CH2 region in the density profiles (Figure 4). 3.4. Micelle Formation. With C36-UA’s reduction in atom interaction sites compared to the AA C36 force field, the time scale for micelle formation can be reached in a more computationally efficient manner. In total, six simulations, M1−M6, were run (Table 2), all above the critical micelle concentration of 1 mM for DPC,53 which yielded a total of twelve micelles. Initially, surfactants formed aggregates which continued to grow (Figure S8 of the Supporting Information), taking up free monomers (Figure 6B) until micelles were formed, which were either stable for the remainder of the simulation or combined again (Figure 6A). All six simulations featured the formation of micelles. The comparison of results with experiments and other simulations will be discussed in Section 4. These results include the shape of micelles, which can be quantified by several properties: aggregation number, radius of gyration (Rg), effective radius (rm), radial distribution functions (RDFs), and the ratio of the maximum to the minimum moment of inertia (Imax/Imin). The aggregation number is the number of surfactants in a micelle, taken at the end of the simulation. The radius of gyration of a micelle is the average distance between the center of mass and every atom of the micelle; no mass weighting is applied. The effective radius is given by54 rm = rh + d − rw

Figure 4. Symmetrized EDPs for DMPC/cholesterol bilayers with C36-UA and C36.4

(2)

rh is the average distance between the center of mass of the micelle and every phosphorus. This was averaged over each time step between the last 50 ns of the simulation. d is the yvalue of the first peak in the O−P RDF of water and phosphorus for a given micelle. rw is the radius of water, 1.4 Å, which represents the relatively immobile shell of hydration around the zwitterionic head groups. Imax/Imin describes the shape of the micelle. A larger Imax/Imin indicates a more elliptical shape, while the lowest limit of unity indicates an ideal sphere. The aforementioned structural characteristics are presented in Table 5. Results from simulations by Cheng et al.55 are included in column 4 to compare the results of the C36-UA simulations here to those of AA C36 simulations; this

of the Supporting Information), the density for C36-UA for the methyl hydrocarbon region (z = 5−12 Å) was higher than that of C36 because of the OPLS CH2 parameters. The deviation was less for DMPC/30% cholesterol. However, consistent with the lower SA/lipid, the head-to-head spacing was larger for C36-UA. The NMR SCDs are compared for the DMPC/cholesterol bilayers in Figure 5. The results were nearly identical for the 30% cholesterol bilayers, and C36-UA had a slightly lower SCD compared to that of C36 for the 10% cholesterol bilayer. 551



Article

Figure 6. (A) Fusion of two micelles from M1. (B) Integration of a lone surfactant to a micelle in M4. (C) Micelles from M3 and M5 at 148.0 ns.

Table 5. Micelle Properties from This Self-Assembly Work and Those of Preformed Micelles by Cheng et al.55,a Imax/Imin rm [Å] Rg [Å] no. of DPC no. of water no. of ions no. of atoms temp (K) no. of micelles [DPC] (mol/L)

M1−M6

M1

M2 and M3

Cheng et al.55

1.32 ± 0.12 16.4 ± 0.4 15.5 ± 0.3 35 ± 18 (14−61) 13 684 none 43 713 300 12 0.30 ± 0.16 (0.1−0.45)

1.23 ± 0.08 18.9 ± 0.2 17.7 ± 0.2 54 6 667 none 22 053 300 1 0.45

1.25 ± 0.02 18.9 ± 0.1 17.8 ± 0.1 53 ± 11 (42−64) 16 667 none 55 130 300 4 0.45

1.22 ± 0.04 21.4 ± 0.1 16.5 ± 0.1 54 11 988 23 K+ and Cl− 39 304 303.15 5 0.25

a For reference, experimental estimates in rm are 19.5−24.553 and a previous MD simulation had an rm of 21−22 Å and an Rg of 17.4 Å.56 M2 and M3 contained only completely equilibrated micelles all greater than 40 DPC surfactants. The reported no. of DPC row has ± standard deviation with the range given in parentheses.

In each simulation, the number of free monomer surfactants approached zero or a small number at long time scales, as represented by the RDF. The formation of a micelle (Figure S9 of the Supporting Information) or the fusion of two micelles (Figure S10 of the Supporting Information) changes the shape of the number density RDF. A number density RDF representing one micelle was roughly sigmoidal (Figure S9). A number density RDF representing two micelles will have points of inflection surrounding a region where the slope decreases almost to a plateau, representing the lack of surfactants between the micelles (Figure S11 of the Supporting Information). Before the fusion (40 ns), M1 had a region of reduced slope (between 20 and 25 Å) and was sigmoidal afterward (Figure S10 of the Supporting Information). M3 approached a sigmoidal shape as the micelle incorporated more of the free monomers (Figure S12 of the Supporting Information). The number density RDF of simulations with multiple micelles (M2 and M4−M6) is more difficult to analyze because the influence of micelles overlap. The length of time required for micelle equilibration is an important characteristic of micelle formation as it can aid researchers in planning and designing simulations.55 An indication that the system is at equilibrium at some time point is that the number density RDF does not change appreciably thereafter (i.e., no major changes in the system structure). To check for equilibrium, it is useful to view the progression of the entire simulation (150 ns) with changes over

comparison will be discussed further in Section 4. Column 1 contains the average for all six simulations. Column 2 contains only results from M1 which has the same number of surfactants as those in the simulations by Cheng et al.55 Column 3 shows the average of only pseudoequilibrated micelles. Unfortunately, not all of the simulations are equally useful because some micelles are not likely in aggregation equilibrium as a bulk solution of DPC molecules. Finite size effects influence the ability to self-aggregate, and simulations M1, M4, and M6 resulted in a single micelle of all DPC molecules available (Table S1 of the Supporting Information). This suggests additional surfactants could be added to the micelle if they were available, but the finite size of the simulation prevented this from occurring. M5 formed micelles of relatively small aggregation numbers of 17 and 20 DPC molecules that, if run longer, may combine into a single micelle. Therefore, only M2 and M3 contained micelles that may be naturally occurring in bulk solution. For our analysis in Table 5, only micelles of aggregation numbers greater than 40 were considered equilibrium micelles, as the experimental aggregation number is 56 ± 5.53 Various RDFs aid analysis of the dynamics of micelle formation. Three atom pairs were used for RDF creation: P−P (headgroup phosphorus), C212−C212 (carbon at the end of the chain, opposite to the headgroup), and C26−C26 (carbon in the middle of the chain). 552



Article

SA/lipid agreed with the experimental estimate (Table 4) that was estimated primarily by the neutron data and (2) DHH agreed reasonably well with the SDP model which gives the strongest signal in the X-ray data. The fine points regarding component densities may be difficult to obtain with scattering densities alone. MD simulations of the C36 and C36-UA bilayers further showed (originally discussed in ref 39) that obtaining certain component distributions with the SDP model (specifically CH) are difficult without constraining fits with the SDP models. There was a consistent discrepancy in the placement of the CH density, i.e., the SDP model suggested a broader distribution shifted away from the bilayer by ∼3 Å compared to that of the MD simulations (Figure 2). The SDP model placed the peak of the CH2 density near 6 Å which is nearly identical to the location of the peak of the CH density from MD. Because the direct comparison with the experimental form factors is the total density, the total chain density for C36 matched the SDP model (Figure S14 of the Supporting Information). It is likely that a similar fit to the SDP model with the CH density closer to the center of the bilayer than that presented by Kučerka et al.39 may be equally valid. More recent SDP analyses have constrained the CH peak to 7.24 Å in agreement with MD simulations using the GROMOS 43A1-S3 force field.42 It is clear from ref 42 and again showed here that either the SDP model alone cannot be used to determine some of the component density locations and experimental probes or MD simulations should aid in constraining the SDP model. Mixing AA and UA representations in the same MD simulation worked well for the DMPC(UA)/cholesterol(AA) test case. The increased CH2 density persisted in these mixtures (Figure 4). However, this did not impact the SCD’s for these mixtures, and the agreement with experiment was comparable (Figure 5). Although we do not expect major issues with the energetic representation of systems using a mixture of the UA-chain and AA CHARMM protein force fields, this test was deemed beyond the scope of this initial project. Micelles successfully self-assembled from an initially random configuration of DPC and water using C36-UA. In the simulations from this work, micelles of a smaller aggregation number (

CHARMM36 United Atom Chain Model for Lipids and Surfactants

Recommend Documents