Ion-Induced PIP2 Clustering with Martini3: Modification of Phosphate-Ion Interactions and Comparison with CHARMM36.

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), K[DMP], and Ca[DMP]2 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.

Spontaneous curvature generation by peptides in asymmetric bilayers.

Peptides and proteins play crucial roles in membrane remodeling by inducing spontaneous curvature. However, extracting spontaneous curvatures from simulations of asymmetric bilayers is challenging because differential stress (i.e., the difference of the leaflet surface tensions) arising from leaflet area strains can vary substantially among initial conditions. This study investigates peptide-induced spontaneous curvature δc 0 p in asymmetric bilayers consisting of a single lipid type and a peptide confined to one leaflet; δc 0 p is calculated from the Helfrich equation using the first moment of the lateral pressure tensor and an alternative expression using the differential stress. It is shown that differential stress introduced during initial system generation is effectively relaxed by equilibrating using P21 periodic boundary conditions, which allows lipids to switch leaflets across cell boundaries and equalize their chemical potentials across leaflets. This procedure leads to robust estimates of δc 0 p for the systems simulated, and is recommended when equality of chemical potentials between the leaflets is a primary consideration.

Flipped C-Terminal Ends of APOA1 Promote ABCA1-Dependent Cholesterol Efflux by Small HDLs.

BACKGROUND: Cholesterol efflux capacity (CEC) predicts cardiovascular disease independently of high-density lipoprotein (HDL) cholesterol levels. Isolated small HDL particles are potent promoters of macrophage CEC by the ABCA1 (ATP-binding cassette transporter A1) pathway, but the underlying mechanisms are unclear. METHODS: We used model system studies of reconstituted HDL and plasma from control and lecithin-cholesterol acyltransferase (LCAT)-deficient subjects to investigate the relationships among the sizes of HDL particles, the structure of APOA1 (apolipoprotein A1) in the different particles, and the CECs of plasma and isolated HDLs. RESULTS: We quantified macrophage and ABCA1 CEC of 4 distinct sizes of reconstituted HDL. CEC increased as particle size decreased. Tandem mass spectrometric analysis of chemically cross-linked peptides and molecular dynamics simulations of APOA1, the major protein of HDL, indicated that the mobility of C-terminus of that protein was markedly higher and flipped off the surface in the smallest particles. To explore the physiological relevance of the model system studies, we isolated HDL from LCAT-deficient subjects, whose small HDLs (like reconstituted HDLs) are discoidal and composed of APOA1, cholesterol, and phospholipid. Despite their very low plasma levels of HDL particles, these subjects had normal CEC. In both the LCAT-deficient subjects and control subjects, the CEC of isolated extra-small HDL (a mixture of extra-small and small HDL by calibrated ion mobility analysis) was 3- to 5-fold greater than that of the larger sizes of isolated HDL. Incubating LCAT-deficient plasma and control plasma with human LCAT converted extra-small and small HDL particles into larger particles, and it markedly inhibited CEC. CONCLUSIONS: We present a mechanism for the enhanced CEC of small HDLs. In smaller particles, the C-termini of the 2 antiparallel molecules of APOA1 are "flipped" off the lipid surface of HDL. This extended conformation allows them to engage with ABCA1. In contrast, the C-termini of larger HDLs are unable to interact productively with ABCA1 because they form a helical bundle that strongly adheres to the lipid on the particle. Enhanced CEC, as seen with the smaller particles, predicts decreased cardiovascular disease risk. Thus, extra-small and small HDLs may be key mediators and indicators of the cardioprotective effects of HDL.

Flipped C-Terminal Ends of APOA1 Promote ABCA1-dependent Cholesterol Efflux by Small HDLs.

Background: Cholesterol efflux capacity (CEC) predicts cardiovascular disease (CVD) independently of HDL cholesterol (HDL-C) levels. Isolated small HDL particles are potent promoters of macrophage CEC by the ABCA1 pathway, but the underlying mechanisms are unclear. Methods: We used model system studies of reconstituted HDL and plasma from control and lecithin-cholesterol acyltransferase (LCAT)-deficient subjects to investigate the relationships among the sizes of HDL particles, the structure of APOA1 in the different particles, and the CECs of plasma and isolated HDLs. Results: We quantified macrophage and ABCA1 CEC of four distinct sizes of reconstituted HDL (r-HDL). CEC increased as particle size decreased. MS/MS analysis of chemically crosslinked peptides and molecular dynamics simulations of APOA1 (HDL's major protein) indicated that the mobility of that protein's C-terminus was markedly higher and flipped off the surface in the smallest particles. To explore the physiological relevance of the model system studies, we isolated HDL from LCAT-deficient subjects, whose small HDLs-like r-HDLs-are discoidal and composed of APOA1, cholesterol, and phospholipid. Despite their very low plasma levels of HDL particles, these subjects had normal CEC. In both the LCAT-deficient subjects and control subjects, the CEC of isolated extra-small HDL (a mixture of extra-small and small HDL by calibrated ion mobility analysis) was 3-5-fold greater than that of the larger sizes of isolated HDL. Incubating LCAT-deficient plasma and control plasma with human LCAT converted extra-small and small HDL particles into larger particles, and it markedly inhibited CEC. Conclusions: We present a mechanism for the enhanced CEC of small HDLs. In smaller particles, the C-termini of the two antiparallel molecules of APOA1 are flipped off the lipid surface of HDL. This extended conformation allows them to engage with ABCA1. In contrast, the C-termini of larger HDLs are unable to interact productively with ABCA1 because they form a helical bundle that strongly adheres to the lipid on the particle. Enhanced CEC, as seen with the smaller particles, predicts decreased CVD risk. Thus, extra-small and small HDLs may be key mediators and indicators of HDL's cardioprotective effects.

Short hydrocarbon stapled ApoC2-mimetic peptides activate lipoprotein lipase and lower plasma triglycerides in mice.

Introduction: Defects in lipolysis can lead to hypertriglyceridemia, which can trigger acute pancreatitis and is also associated with cardiovascular disease. Decreasing plasma triglycerides (TGs) by activating lipoprotein lipase (LPL) with ApoC2 mimetic peptides is a new treatment strategy for hypertriglyceridemia. We recently described a dual ApoC2 mimetic/ApoC3 antagonist peptide called D6PV that effectively lowered TG in several mouse models but has limitations in terms of drug development. The aim of this study was to create the next generation of ApoC2 mimetic peptides. Methods: We employed hydrocarbon staples, as well as select amino acid substitutions, to make short single helical mimetic peptides based on the last helix of ApoC2. Peptides were first tested for their ability to activate LPL and then in hypertriglyceridemia mouse models. All-atom simulations of peptides were performed in a lipid-trilayer model of TG-rich lipoproteins to discern their possible mechanism of action. Results: We designed a single stapled peptide called SP1 (21 residues), and a double stapled (stitched) peptide called SP2 (21 residues) and its N-terminal acylated analogue, SP2a. The hydrocarbon staples increased the amphipathicity of the peptides and their ability to bind lipids without interfering with LPL activation. Indeed, from all-atom simulations, the conformations of SP1 and SP2a are restrained by the staples and maintains the proper orientation of the LPL activation motif, while still allowing their deeper insertion into the lipid-trilayer model. Intraperitoneal injection of stapled peptides (1-5 umoles/kg) into ApoC2-hypomorphic mice or human ApoC3-transgenic resulted in an 80%-90% reduction in plasma TG within 3 h, similar to the much longer D6PV peptide (41 residues). Other modifications (replacement L-Glu20, L-Glu21 with their D-isomers, N-methylation of Gly19, Met2NorLeu and Ala1alpha-methylAla substitutions, N-terminal octanoylation) were introduced into the SP2a peptide. These changes made SP2a highly resistant to proteolysis against trypsin, pepsin, and Proteinase K, while maintaining similar efficacy in lowering plasma TG in mice. Conclusion: We describe a new generation of ApoC2 mimetic peptides based on hydron carbon stapling that are at least equally potent to earlier peptides but are much shorter and resistant to proteolysis and could be further developed into a new therapy for hypertriglyceridemia.

Unsaturated Lipids Facilitate Partitioning of Transmembrane Peptides into the Liquid Ordered Phase.

The affinity of single-pass transmembrane (TM) proteins for ordered membrane phases has been reported to depend on their lipidation, TM length, and lipid accessible surface area. In this work, the raft affinities of the TM domain of the linker for activation of T cells and its depalmitoylated variant are assessed using free energy simulations in a binary bilayer system composed of two laterally patched bilayers of ternary liquid ordered (Lo) and liquid disordered (Ld) phases. These phases are modeled by distinct compositions of distearoylphosphatidylcholine, palmitoyloleoylphosphatidylcholine (POPC), and cholesterol, and the simulations were carried out for 4.5 µs/window. Both peptides are shown to preferentially partition into the Ld phase in agreement with model membrane experiments and previous simulations on ternary lipid mixtures but not with measurements on giant plasma membrane vesicles where the Lo is slightly preferred. However, the 500 ns average relaxation time of lipid rearrangement around the peptide precluded a quantitative analysis of free energy differences arising from peptide palmitoylation and two different lipid compositions. When in the Lo phase, peptides reside in regions rich in POPC and interact preferentially with its unsaturated tail. Hence, the detailed substructure of the Lo phase is an important modulator of peptide partitioning, in addition to the inherent properties of the peptide.

Drude Polarizable Lipid Force Field with Explicit Treatment of Long-Range Dispersion: Parametrization and Validation for Saturated and Monounsaturated Zwitterionic Lipids.

Accurate empirical force fields of lipid molecules are a critical component of molecular dynamics simulation studies aimed at investigating properties of monolayers, bilayers, micelles, vesicles, and liposomes, as well as heterogeneous systems, such as protein-membrane complexes, bacterial cell walls, and more. While the majority of lipid force field-based simulations have been performed using pairwise-additive nonpolarizable models, advances have been made in the development of the polarizable force field based on the classical Drude oscillator model. In the present study, we undertake further optimization of the Drude lipid force field, termed Drude2023, including improved treatment of the phosphate and glycerol linker region of PC and PE headgroups, additional optimization of the alkene group in monounsaturated lipids, and inclusion of long-range Lennard-Jones interactions using the particle-mesh Ewald method. Initial optimization targeted quantum mechanical (QM) data on small model compounds representative of the linker region. Subsequent optimization targeted QM data on larger model compounds, experimental data, and dihedral potentials of mean force from the CHARMM36 additive lipid force field using a parameter reweighting protocol. The use of both experimental and QM target data during the reweighting protocol is shown to produce physically reasonable parameters that reproduce a collection of experimental observables. Target data for optimization included surface area/lipid for DPPC, DSPC, DMPC, and DLPC bilayers and nuclear magnetic resonance (NMR) order parameters for DPPC bilayers. Validation data include prediction of membrane thickness, scattering form factors, electrostatic potential profiles, compressibility moduli, surface area per lipid, water permeability, NMR T1 relaxation times, diffusion constants, and monolayer surface tensions for a variety of saturated and unsaturated lipid mono- and bilayers. Overall, the agreement with experimental data is quite good, though the results are less satisfactory for the NMR T1 relaxation times for carbons near the ester groups. Notable improvements compared to the additive C36 force field were obtained for membrane dipole potentials, lipid diffusion coefficients, and water permeability with the exception of monounsaturated lipid bilayers. It is anticipated that the optimized polarizable Drude2023 force field will help generate more accurate molecular simulations of pure bilayers and heterogeneous systems containing membranes, advancing our understanding of the role of electronic polarization in these systems.

Surface viscosities of lipid bilayers determined from equilibrium molecular dynamics simulations.

Lipid membrane viscosity is critical to biological function. Bacterial cells grown in different environments alter their lipid composition in order to maintain a specific viscosity, and membrane viscosity has been linked to the rate of cellular respiration. To understand the factors that determine the viscosity of a membrane, we ran equilibrium all-atom simulations of single component lipid bilayers and calculated their viscosities. The viscosity was calculated via a Green-Kubo relation, with the stress-tensor autocorrelation function modeled by a stretched exponential function. By simulating a series of lipids at different temperatures, we establish the dependence of viscosity on several aspects of lipid chemistry, including hydrocarbon chain length, unsaturation, and backbone structure. Sphingomyelin is found to have a remarkably high viscosity, roughly 20 times that of DPPC. Furthermore, we find that inclusion of the entire range of the dispersion interaction increases viscosity by up to 140%. The simulated viscosities are similar to experimental values obtained from the rotational dynamics of small chromophores and from the diffusion of integral membrane proteins but significantly lower than recent measurements based on the deformation of giant vesicles.

Initiation and evolution of pores formed by influenza fusion peptides probed by lysolipid inclusion.

The fusion peptide (FP) domain is necessary for the fusogenic activity of spike proteins in a variety of enveloped viruses, allowing the virus to infect the host cell, and is the only part of the protein that interacts directly with the target membrane lipid tails during fusion. There are consistent findings of poration by this domain in experimental model membrane systems, and, in certain conditions, the isolated FPs can generate pores. Here, we use molecular dynamics simulations to investigate the specifics of how these FP-induced pores form in membranes with different compositions of lysolipid and POPC. The simulations show that pores form spontaneously at high lysolipid concentrations via hybrid intermediates, where FP aggregates in the cis leaflet tilt to form a funnel-like structure that spans the leaflet and locally reduces the hydrophobic thickness that must be traversed by water to form a pore. By restraining a single FP within an FP aggregate to this tilted conformation, pores can be formed in lower-lysolipid-content membranes, including pure POPC, on the 100-ns timescale, much more rapidly than in unbiased simulations in bilayers with the same composition. The pore formation pathway is similar to the spontaneous formation in high lysolipid concentrations. Depending on the membrane composition, the pores can be metastable (as seen in POPC) or lead to membrane rupture.

Planar aggregation of the influenza viral fusion peptide alters membrane structure and hydration, promoting poration.

To infect, enveloped viruses employ spike protein, spearheaded by its amphipathic fusion peptide (FP), that upon activation extends out from the viral surface to embed into the target cellular membrane. Here we report that synthesized influenza virus FPs are membrane active, generating pores in giant unilamellar vesicles (GUV), and thus potentially explain both influenza virus' hemolytic activity and the liposome poration seen in cryo-electron tomography. Experimentally, FPs are heterogeneously distributed on the GUV at the time of poration. Consistent with this heterogeneous distribution, molecular dynamics (MD) simulations of asymmetric bilayers with different numbers of FPs in one leaflet show FP aggregation. At the center of FP aggregates, a profound change in the membrane structure results in thinning, higher water permeability, and curvature. Ultimately, a hybrid bilayer nanodomain forms with one lipidic leaflet and one peptidic leaflet. Membrane elastic theory predicts a reduced barrier to water pore formation when even a dimer of FPs thins the membrane as above, and the FPs of that dimer tilt, to continue the leaflet bending initiated by the hydrophobic mismatch between the FP dimer and the surrounding lipid.

Finite-Size Effects in Simulations of Peptide/Lipid Assembly.

Molecular dynamics simulations are an attractive tool for understanding lipid/peptide self-assembly but can be plagued by inaccuracies when the system sizes are too small. The general guidance from self-assembly simulations of homogeneous micelles is that the total number of surfactants should be three to five times greater than the equilibrium aggregate number of surfactants per micelle. Herein, the heuristic is tested on the more complicated self-assembly of lipids and amphipathic peptides using the Cooke and Martini 3 coarse-grained models. Cooke model simulations with 50 to 1000 lipids and no peptide are dominated by finite-size effects, with usually one aggregate (micelle or nanodisc) containing most of the lipids forming at each system size. Approximately 200 systems of different peptide/lipid (P/L) ratios and sizes of up to 1000 lipids yield a "finite-size phase diagram" for peptide driven self-assembly, including a coexistence region of micelles and discs. Insights from the Cooke model are applied to the assembly of dimyristoylphosphatidylcholine and the ELK-neutral peptide using the Martini 3 model. Systems of 150, 450, and 900 lipids with P/L = 1/6.25 form mixtures of lipid-rich discs that agree in size with experiment and peptide-rich micelles. Only the 150-lipid system shows finite-size effects, which arise from the long-tailed distribution of aggregate sizes. The general rule of three to five times the equilibrium aggregate size remains a practical heuristic for the Cooke and Martini 3 systems investigated here.

Design principles of PI(4,5)P2 clustering under protein-free conditions: Specific cation effects and calcium-potassium synergy.

Phosphatidylinositol 4,5-bisphosphate (PIP2) clustering is a key component in cell signaling, yet little is known about the atomic-level features of this phenomenon. Network-theoretic analysis of multimicrosecond atomistic simulations of PIP2 containing asymmetric bilayers under protein-free conditions, presented here, reveals how design principles of PIP2 clustering are determined by the specific cation effects. Ca2+ generates large clusters (6% are pentamer or larger) by adding existing PIP2 dimers formed by strong O‒Ca2+‒O bridging interactions of unprotonated P4/P5 phosphates. In contrast, monovalent cations (Na+ and K+) form smaller and less-stable clusters by preferentially adding PIP2 monomers. Despite having the same net charge, the affinity to P4/P5 is higher for Na+, while affinity toward glycerol P1 is higher for K+. Consequently, a mixture of K+ and Ca2+ (as would be produced by Ca2+ influx) synergistically yields larger and more stable clusters than Ca2+ alone due to the different binding preferences of these cations.

Development of CHARMM Additive Potential Energy Parameters for α-Methyl Amino Acids.

Potential energy parameters for α-methyl amino acids were generated with ab initio calculations on α-methyl-N-acetylalanyl-N'-methylamide (the α-methyl "alanine dipeptide") which served as an input to a grid-based correction to the backbone torsional potential (known as CMAP) consistent with the CHARMM36m additive protein force field. The new parameters were validated by comparison with experimentally determined helicities of the 22 residue C-terminal peptide (H10) from apolipoprotein A1 and five α-methylated variants in water and 0.3:0.7 trifluoroethanol (TFE)/water. Conventional molecular dynamics simulation totaling 30 µs for each peptide is in overall good agreement with the experiment, including the increased helicity in 30% TFE. An additional 500 ns of simulation using two-dimensional dihedral biasing (bpCMAP) replica exchange reduced left-handed conformations, increased right-handed helices, and thereby mostly decreased agreement with the experiment. Analysis of side chain-side chain salt bridges suggests that the overestimation of the helical content may be, in part, due to such interactions. The increased helicity of the peptides in 30% TFE arises from decreased hydrogen bonding of the backbone atoms to water and a concomitant increase in intramolecular backbone hydrogen bonds.

Developing initial conditions for simulations of asymmetric membranes: a practical recommendation.

It has been proposed that the surface tension difference between leaflets (or differential stress) in asymmetric bilayers is generally nonvanishing. This implies that there is no unique approach to generate initial conditions for simulations of asymmetric bilayers in the absence of experimentally derived constraints. Current generation methods include individual area per lipid (APL) based, leaflet surface area (SA) matching, and zero leaflet tension based (0-DS). This work adds a bilayer-based approach that aims for achieving partial chemical equilibrium by interleaflet switching of selected lipids via P21 periodic boundary conditions. Based on a recently proposed theoretical framework, we obtained expressions for tensions in asymmetric bilayers from both the bending and area strains. We also developed a quantitative measure for the energetic penalty from the differential stress. The impacts of APL-, SA-, and 0-DS-based approaches on mechanical properties are assessed for two different asymmetric bilayers. The lateral pressure profile and its moments differ significantly for each method, whereas the area compressibility modulus is relatively insensitive. Application of P21 periodic boundary conditions (APL/P21, SA/P21, and 0-DS/P21) results in better agreement in mechanical properties between asymmetric bilayers generated by APL-, SA-, and 0-DS-based approaches, in which changes are the smallest for bilayers from the SA-based method. The estimated differential stress from the theory shows good agreement with that from the simulations. These simulation results and the good agreement between the predicted and observed differential stress further support the theoretical framework in which bilayer mechanical properties are outcomes of the interplay between intrinsic bending and asymmetric lipid packing. Based on the simulation results and theoretical predictions, the SA/P21-based, or at least the SA-based (when the differential stress is small), approach is recommended as a practical method for developing initial conditions for asymmetric bilayer simulations.

Location and Conformational Ensemble of Menaquinone and Menaquinol, and Protein-Lipid Modulations in Archaeal Membranes.

Halobacteria, a type of archaea in high salt environments, have phytanyl ether phospholipid membranes containing up to 50% menaquinone. It is not understood why a high concentration of menaquinone is required and how it influences membrane properties. In this study, menaquinone-8 headgroup and torsion parameters of isoprenoid tail are optimized in the CHARMM36 force field. Molecular dynamics simulations of archaeal bilayers containing 0 to 50% menaquinone characterize the distribution of menaquinone-8 and menaquinol-8, as well as their effects on mechanical properties and permeability. Menaquinone-8 segregates to the membrane midplane above concentrations of 10%, favoring an extended conformation in a fluid state. Menaquinone-8 increases the bilayer thickness but does not significantly alter the area compressibility modulus and lipid chain ordering. Counterintuitively, menaquinone-8 increases water permeability because it lowers the free energy barrier in the midplane. The thickness increase due to menaquinone-8 may help halobacteria ameliorate hyper-osmotic pressure by increasing the membrane bending constant. Simulations of the archaeal membranes with archaerhodopsin-3 show that the local membrane surface adjusts to accommodate the thick membranes. Overall, this study delineates the biophysical landscape of 50% menaquinone in the archaeal bilayer, demonstrates the mixing of menaquinone and menaquinol, and provides atomistic details about menaquinone configurations.

Functional Group Distributions, Partition Coefficients, and Resistance Factors in Lipid Bilayers Using Site Identification by Ligand Competitive Saturation.

Small molecules such as metabolites and drugs must pass through the membrane of the cell, a barrier primarily comprising phospholipid bilayers and embedded proteins. To better understand the process of passive diffusion, knowledge of the ability of various functional groups to partition across bilayers and the associated energetics would be of utility. In the present study, the site identification by ligand competitive saturation (SILCS) methodology has been applied to sample the distributions of a diverse set of chemical solutes representing the functional groups of small molecules across phospholipid bilayers composed of 0.9:0.1 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine/cholesterol and a mixture of 0.52:0.18:0.3 1,2-dioleoyl-sn-glycero-3-phospho-l-serine/1,2-dioleoyl-sn-glycero-3-phosphocholine/cholesterol used in parallel artificial membrane permeability assay experiments. A combination of oscillating chemical potential grand canonical Monte Carlo and molecular dynamics in the SILCS simulations was applied to achieve solute sampling through the bilayers and surrounding aqueous environment from which the distribution of solutes and the functional groups they represent were obtained. Results show differential distribution of aliphatic versus aromatic groups with the former having increased sampling in the center of the bilayers versus in the region of the glycerol linker for the latter. Variations in the distribution of different polar groups are evident, with large differences between negative acetate and positive methylammonium with accumulation of the polar-neutral and acetate solutes above the bilayer head groups. Conversion of the distributions to absolute free energies allows for a detailed understanding of energetics of functional groups in different regions of the bilayers and for calculation of absolute free-energy profiles of multifunctional drug-like molecules across the bilayers from which partition coefficients and resistance factors suitable for insertion into the homogenous solubility-diffusion equation for calculation of permeability were obtained. Comparisons of the calculated bilayer/solution partition coefficients with 1-octanol/water experimental data for both drug-like molecules and the solutes show overall good agreement, validating the calculated distributions and associated absolute free-energy profiles.

CHARMM36 Lipid Force Field with Explicit Treatment of Long-Range Dispersion: Parametrization and Validation for Phosphatidylethanolamine, Phosphatidylglycerol, and Ether Lipids.

Long-range Lennard-Jones (LJ) interactions have been incorporated into the CHARMM36 (C36) lipid force field (FF) using the LJ particle-mesh Ewald (LJ-PME) method in order to remove the inconsistency of bilayer and monolayer properties arising from the exclusion of long-range dispersion [Yu, Y.; Semi-automated Optimization of the CHARMM36 Lipid Force Field to Include Explicit Treatment of Long-Range Dispersion. J. Chem. Theory Comput. 2021, 10.1021/acs.jctc.0c01326. (preceding article in this issue)]. The new FF is denoted C36/LJ-PME. While the first optimization was based on three phosphatidylcholines (PCs), this work extends the validation and parametrization to more lipids including PC, phosphatidylethanolamine (PE), phosphatidylglycerol (PG), and ether lipids. The agreement with experimental structure data is excellent for PC, PE, and ether lipids. C36/LJ-PME also compares favorably with scattering data of PG bilayers but less so with NMR deuterium order parameters of 1,2-dimyristoyl-sn-glycero-3-phospho-(1'-rac-glycerol) (DMPG) at 303.15 K, indicating a need for future optimization regarding PG-specific parameters. Frequency dependence of NMR T1 spin-lattice relaxation times is well-described by C36/LJ-PME, and the overall agreement with experiment is comparable to C36. Lipid diffusion is slower than C36 due to the added long-range dispersion causing a higher viscosity, although it is still too fast compared to experiment after correction for periodic boundary conditions. When using a 10 Å real-space cutoff, the simulation speed of C36/LJ-PME is roughly equal to C36. While more lipids will be incorporated into the FF in the future, C36/LJ-PME can be readily used for common lipids and extends the capability of the CHARMM FF by supporting monolayers and eliminating the cutoff dependence.

Semi-automated Optimization of the CHARMM36 Lipid Force Field to Include Explicit Treatment of Long-Range Dispersion.

The development of the CHARMM lipid force field (FF) can be traced back to the early 1990s with its current version denoted CHARMM36 (C36). The parametrization of C36 utilized high-level quantum mechanical data and free energy calculations of model compounds before parameters were manually adjusted to yield agreement with experimental properties of lipid bilayers. While such manual fine-tuning of FF parameters is based on intuition and trial-and-error, automated methods can identify beneficial modifications of the parameters via their sensitivities and thereby guide the optimization process. This work introduces a semi-automated approach to reparametrize the CHARMM lipid FF with consistent inclusion of long-range dispersion through the Lennard-Jones particle-mesh Ewald (LJ-PME) approach. The optimization method is based on thermodynamic reweighting with regularization with respect to the C36 set. Two independent optimizations with different topology restrictions are presented. Targets of the optimizations are primarily liquid crystalline phase properties of lipid bilayers and the compression isotherm of monolayers. Pair correlation functions between water and lipid functional groups in aqueous solution are also included to address headgroup hydration. While the physics of the reweighting strategy itself is well-understood, applying it to heterogeneous, complex anisotropic systems poses additional challenges. These were overcome through careful selection of target properties and reweighting settings allowing for the successful incorporation of the explicit treatment of long-range dispersion, and we denote the newly optimized lipid force field as C36/LJ-PME. The current implementation of the optimization protocol will facilitate the future development of the CHARMM and related lipid force fields.

Membrane permeability of small molecules from unbiased molecular dynamics simulations.

Permeation of many small molecules through lipid bilayers can be directly observed in molecular dynamics simulations on the nano- and microsecond timescale. While unbiased simulations provide an unobstructed view of the permeation process, their feasibility for computing permeability coefficients depends on various factors that differ for each permeant. The present work studies three small molecules for which unbiased simulations of permeation are feasible within less than a microsecond, one hydrophobic (oxygen), one hydrophilic (water), and one amphiphilic (ethanol). Permeabilities are computed using two approaches: counting methods and a maximum-likelihood estimation for the inhomogeneous solubility diffusion (ISD) model. Counting methods yield nearly model-free estimates of the permeability for all three permeants. While the ISD-based approach is reasonable for oxygen, it lacks precision for water due to insufficient sampling and results in misleading estimates for ethanol due to invalid model assumptions. It is also demonstrated that simulations using a Langevin thermostat with collision frequencies of 1/ps and 5/ps yield oxygen permeabilities and diffusion constants that are lower than those using Nosé-Hoover by statistically significant margins. In contrast, permeabilities from trajectories generated with Nosé-Hoover and the microcanonical ensemble do not show statistically significant differences. As molecular simulations become more affordable and accurate, calculation of permeability for an expanding range of molecules will be feasible using unbiased simulations. The present work summarizes theoretical underpinnings, identifies pitfalls, and develops best practices for such simulations.