This Review illustrates the evaluation of permeability of lipid membranes from molecular dynamics (MD) simulation primarily using water and oxygen as examples. Membrane entrance, translocation, and ...exit of these simple permeants (one hydrophilic and one hydrophobic) can be simulated by conventional MD, and permeabilities can be evaluated directly by Fick’s First Law, transition rates, and a global Bayesian analysis of the inhomogeneous solubility-diffusion model. The assorted results, many of which are applicable to simulations of nonbiological membranes, highlight the limitations of the homogeneous solubility diffusion model; support the utility of inhomogeneous solubility diffusion and compartmental models; underscore the need for comparison with experiment for both simple solvent systems (such as water/hexadecane) and well-characterized membranes; and demonstrate the need for microsecond simulations for even simple permeants like water and oxygen. Undulations, subdiffusion, fractional viscosity dependence, periodic boundary conditions, and recent developments in the field are also discussed. Last, while enhanced sampling methods and increasingly sophisticated treatments of diffusion add substantially to the repertoire of simulation-based approaches, they do not address directly the critical need for force fields with polarizability and multipoles, and constant pH methods.
The C36 CHARMM lipid force field has been extended to include sphingolipids, via a combination of high-level quantum mechanical calculations on small molecule fragments, and validation by extensive ...molecular dynamics simulations on N-palmitoyl and N-stearoyl sphingomyelin. NMR data on these two molecules from several studies in bilayers and micelles played a strong role in the development and testing of the force field parameters. Most previous force fields for sphingomyelins were developed before the availability of the detailed NMR data and relied on x-ray diffraction of bilayers alone for the validation; these are shown to be too dense in the bilayer plane based on published chain order parameter data from simulations and experiments. The present simulations reveal O-H:::O-P intralipid hydrogen bonding occurs 99% of the time, and interlipid N-H:::O=C (26-29%, depending on the lipid) and N-H:::O-H (17–19%). The interlipid hydrogen bonds are long lived, showing decay times of 50 ns, and forming strings of lipids, and leading to reorientational correlation time of nearly 100 ns. The spontaneous radius of curvature for pure N-palmitoyl sphingomyelin bilayers is estimated to be 43–100 Å, depending on the assumptions made in assigning a bending constant; this unusual positive curvature for a two-tailed neutral lipid is likely associated with hydrogen bond networks involving the NH of the sphingosine group.
•Mechanical properties of 12 different bilayers are evaluated from simulation.•The CHARMM 36 force field yields excellent area and compressibility.•The relation of surface area and NMR order ...parameters does not fit simple models.•Bending constants equal those from flicker experiments, but not others.•Spontaneous curvatures of leaflets in bilayers are 30% lower than in hexagonal phases.
Lipid areas (Aℓ), bilayer area compressibilities (KA), bilayer bending constants (KC), and monolayer spontaneous curvatures (c0) from simulations using the CHARMM36 force field are reported for 12 representative homogenous lipid bilayers. Aℓ (or their surrogate, the average deuterium order parameter in the “plateau region” of the chain) agree very well with experiment, as do the KA. Simulated KC are in near quantitative agreement with vesicle flicker experiments, but are somewhat larger than KC from X-ray, pipette aspiration, and neutron spin echo for saturated lipids. Spontaneous curvatures of bilayer leaflets from the simulations are approximately 30% smaller than experimental values of monolayers in the inverse hexagonal phase.
A significant modification to the additive all-atom CHARMM lipid force field (FF) is developed and applied to phospholipid bilayers with both choline and ethanolamine containing head groups and with ...both saturated and unsaturated aliphatic chains. Motivated by the current CHARMM lipid FF (C27 and C27r) systematically yielding values of the surface area per lipid that are smaller than experimental estimates and gel-like structures of bilayers well above the gel transition temperature, selected torsional, Lennard-Jones and partial atomic charge parameters were modified by targeting both quantum mechanical (QM) and experimental data. QM calculations ranging from high-level ab initio calculations on small molecules to semiempirical QM studies on a 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC) bilayer in combination with experimental thermodynamic data were used as target data for parameter optimization. These changes were tested with simulations of pure bilayers at high hydration of the following six lipids: DPPC, 1,2-dimyristoyl-sn-phosphatidylcholine (DMPC), 1,2-dilauroyl-sn-phosphatidylcholine (DLPC), 1-palmitoyl-2-oleoyl-sn-phosphatidylcholine (POPC), 1,2-dioleoyl-sn-phosphatidylcholine (DOPC), and 1-palmitoyl-2-oleoyl-sn-phosphatidylethanolamine (POPE); simulations of a low hydration DOPC bilayer were also performed. Agreement with experimental surface area is on average within 2%, and the density profiles agree well with neutron and X-ray diffraction experiments. NMR deuterium order parameters (S CD) are well predicted with the new FF, including proper splitting of the S CD for the aliphatic carbon adjacent to the carbonyl for DPPC, POPE, and POPC bilayers. The area compressibility modulus and frequency dependence of 13C NMR relaxation rates of DPPC and the water distribution of low hydration DOPC bilayers also agree well with experiment. Accordingly, the presented lipid FF, referred to as C36, allows for molecular dynamics simulations to be run in the tensionless ensemble (NPT), and is anticipated to be of utility for simulations of pure lipid systems as well as heterogeneous systems including membrane proteins.
Overbinding of ions to lipid head groups is a potentially serious artifact in simulations of charged lipid bilayers. In this study, the Lennard-Jones radii in the CHARMM force field for interactions ...of Na+ and lipid oxygen atoms of carboxyl, phosphate, and ester groups were revised to match osmotic pressure data on sodium acetate and electrophoresis data on palmitoyloleoyl phosphatidylcholine (POPC) vesicles. The new parameters were then validated by successfully reproducing previously published experimental NMR deuterium order parameters for dimyristoyl phosphatidylglycerol (DMPG) and newly obtained values for palmitoyloleoyl phosphatidylserine (POPS). Although the increases in Lennard-Jones diameters are only 0.02–0.12 Å, they are sufficient to reduce Na+ binding, and thereby increase surface areas per lipid by 5–10% compared with the unmodified parameters.
The functional significance of ordered nanodomains (or rafts) in cholesterol rich eukaryotic cell membranes has only begun to be explored. This study exploits the correspondence of cellular rafts and ...liquid ordered (L
) phases of three-component lipid bilayers to examine permeability. Molecular dynamics simulations of L
phase dipalmitoylphosphatidylcholine (DPPC), dioleoylphosphatidylcholine (DOPC), and cholesterol show that oxygen and water transit a leaflet through the DOPC and cholesterol rich boundaries of hexagonally packed DPPC microdomains, freely diffuse along the bilayer midplane, and escape the membrane along the boundary regions. Electron paramagnetic resonance experiments provide critical validation: the measured ratio of oxygen concentrations near the midplanes of liquid disordered (L
) and L
bilayers of DPPC/DOPC/cholesterol is 1.75 ± 0.35, in very good agreement with 1.3 ± 0.3 obtained from simulation. The results show how cellular rafts can be structurally rigid signaling platforms while remaining nearly as permeable to small molecules as the L
phase.
We present an extension of the CHARMM hexopyranose monosaccharide additive all-atom force field to enable modeling of glycosidic-linked hexopyranose polysaccharides. The new force field parameters ...encompass 1→1, 1→2, 1→3, 1→4, and 1→6 hexopyranose glycosidic linkages, as well as O-methylation at the C1 anomeric carbon, and are developed to be consistent with the CHARMM all-atom biomolecular force fields for proteins, nucleic acids, and lipids. The parameters are developed in a hierarchical fashion using model compounds containing the key atoms in the full carbohydrates, in particular O-methyl-tetrahydropyran and glycosidic-linked dimers consisting of two molecules of tetrahyropyran or one molecule of tetrahydropyran and one of cyclohexane. Target data for parameter optimization include full two-dimensional energy surfaces defined by the Φ/Ψ glycosidic dihedral angles in the disaccharide analogs, as determined by quantum mechanical MP2/cc-pVTZ single point energies on MP2/6-31G(d) optimized structures (MP2/cc-pVTZ//MP2/6-31G(d)). In order to achieve balanced, transferable dihedral parameters for the Φ/Ψ glycosidic dihedral angles, surfaces for all possible chiralities at the ring carbon atoms involved in the glycosidic linkages are considered, resulting in over 5 000 MP2/cc-pVTZ//MP2/6-31G(d) conformational energies. Also included as target data are vibrational frequencies, pair interaction energies and distances with water molecules, and intramolecular geometries including distortion of the glycosidic valence angle as a function of the glycosidic dihedral angles. The model compound optimized force field parameters are validated on full disaccharides through the comparison of molecular dynamics results to available experimental data. Good agreement is achieved with experiment for a variety of properties including crystal cell parameters and intramolecular geometries, aqueous densities, and aqueous NMR coupling constants associated with the glycosidic linkage. The newly developed parameters allow for the modeling of linear, branched, and cyclic hexopyranose glycosides both alone and in heterogeneous systems including proteins, nucleic acids, and/or lipids when combined with existing CHARMM biomolecular force fields.
A revision (C35r) to the CHARMM ether force field is shown to reproduce experimentally observed conformational populations of dimethoxyethane. Molecular dynamics simulations of 9, 18, 27, and 36-mers ...of polyethylene oxide (PEO) and 27-mers of polyethylene glycol (PEG) in water based on C35r yield a persistence length
λ
=
3.7
Å, in quantitative agreement with experimentally obtained values of 3.7
Å for PEO and 3.8
Å for PEG; agreement with experimental values for hydrodynamic radii of comparably sized PEG is also excellent. The exponent
υ relating the radius of gyration and molecular weight
(
R
g
∝
M
w
υ
)
of PEO from the simulations equals 0.515
±
0.023, consistent with experimental observations that low molecular weight PEG behaves as an ideal chain. The shape anisotropy of hydrated PEO is 2.59:1.44:1.00. The dimension of the middle length for each of the polymers nearly equals the hydrodynamic radius
R
h obtained from diffusion measurements in solution. This explains the correspondence of
R
h and
R
p, the pore radius of membrane channels: a polymer such as PEG diffuses with its long axis parallel to the membrane channel, and passes through the channel without substantial distortion.
Phosphatidylinositol 4,5-bisphosphate (PIP2) is a critical lipid for cellular signaling. The specific phosphorylation of the inositol ring controls protein binding as well as clustering behavior. Two ...popular models to describe ion-mediated clustering of PIP2 are Martini3 (M3) and CHARMM36 (C36). Molecular dynamics simulations of PIP2-containing bilayers in solutions of potassium chloride, sodium chloride, and calcium chloride, and at two different resolutions are performed to understand the aggregation and the model parameters that drive it. The average M3 clusters of PIP2 in bilayers of 1-palmitoyl-2-oleoyl-glycero-3-phosphocholine and PIP2 bilayers in the presence of K+, Na+, or Ca2+ contained 2.2, 2.6, and 6.4 times more PIP2 than C36 clusters, respectively. Indeed, the Ca2+-containing systems often formed a single large aggregate. Reparametrization of the M3 ion–phosphate Lennard-Jones interaction energies to reproduce experimental osmotic pressure of sodium dimethyl phosphate (DMP), KDMP, and CaDMP2 solutions, the same experimental target as C36, yielded comparably sized PIP2 clusters for the two models. Furthermore, C36 and the modified M3 predict similar saturation of the phosphate groups with increasing Ca2+, although the coarse-grained model does not capture the cooperativity between K+ and Ca2+. This characterization of the M3 behavior in the presence of monovalent and divalent ions lays a foundation to study cation/protein/PIP2 clustering.
PDF
