Coarse-grained Lipid Model Research Articles

Elasticity of bilayer lipid membranes from their density correlation function

ABSTRACT We study the density correlation function (DCF) of DPPC lipid bilayers. We compare Molecular Dynamics (MD) results with theoretical predictions obtained with a mesoscopic description, in terms of the lipid membrane elasticity. One key objective of our work is the quantification of the lipid membrane elasticity directly from the DCF, both for the membrane undulations and local membrane thickness. Our method does not require the definition of instantaneous surfaces or internal variables defining lipid orientations. Building on our previous work, here we focus on the intralayer correlations, i.e. the DCF of lipids residing on the same monolayer, by tracking only the position of the phosphorus atoms in a lipid head group. We demonstrate the relevance of the intralayer two-dimensional (2D) correlations to the total DCF. We further show that all-atom (AA) and coarse grained (CG) lipid forcefields, feature distintively different DCFs. The CG forcefield predicts results in good agreement with the mesoscopic predictions, for the entire wavevector range; the AA forcefield (CHARMM36) predict strong peristaltic fluctuations at long wavevectors q ≳ 0.8 nm − 1 , which are absent in the CG lipid model (MARTINI).

Backmapping with Mapping and Isomeric Information.

I present a powerful and flexible backmapping tool named Multiscale Simulation Tool (mstool) that converts a coarse-grained (CG) system into all-atom (AA) resolution and only requires AA to CG mapping and isomeric information (cis/trans/dihedral/chiral). The backmapping procedure includes two simple steps: (a) AA atoms are randomly placed near the corresponding CG beads according to the provided mapping scheme. (b) Energy minimization is performed with two modifications in the AA force field (FF). First, nonbonded interactions are replaced with cosine functions to ensure the numerical stability. Second, additional torsions are imposed to maintain the molecules' isomeric properties. To test the simplicity and robustness of the tool, I backmapped multiple membrane and protein CG structures into AA resolution, including a four-bead CG lipid model (resolution increased by a factor of 34) without using intermediate resolution. The tool is freely available at github.com/ksy141/mstool.

Automatic multi-objective optimization of coarse-grained lipid force fields using SwarmCG.

The development of coarse-grained (CG) molecular models typically requires a time-consuming iterative tuning of parameters in order to have the approximated CG models behave correctly and consistently with, e.g., available higher-resolution simulation data and/or experimental observables. Automatic data-driven approaches are increasingly used to develop accurate models for molecular dynamics simulations. However, the parameters obtained via such automatic methods often make use of specifically designed interaction potentials and are typically poorly transferable to molecular systems or conditions other than those used for training them. Using a multi-objective approach in combination with an automatic optimization engine (SwarmCG), here, we show that it is possible to optimize CG models that are also transferable, obtaining optimized CG force fields (FFs). As a proof of concept, here, we use lipids for which we can avail reference experimental data (area per lipid and bilayer thickness) and reliable atomistic simulations to guide the optimization. Once the resolution of the CG models (mapping) is set as an input, SwarmCG optimizes the parameters of the CG lipid models iteratively and simultaneously against higher-resolution simulations (bottom-up) and experimental data (top-down references). Including different types of lipid bilayers in the training set in a parallel optimization guarantees the transferability of the optimized lipid FF parameters. We demonstrate that SwarmCG can reach satisfactory agreement with experimental data for different resolution CG FFs. We also obtain stimulating insights into the precision-resolution balance of the FFs. The approach is general and can be effectively used to develop new FFs and to improve the existing ones.

Modeling lipid-protein interactions for coarse-grained lipid and Cα protein models.

Biological membranes that play major roles in diverse functions are composed of numerous lipids and proteins, making them an important target for coarse-grained (CG) molecular dynamics (MD) simulations. Recently, we have developed the CG implicit solvent lipid force field (iSoLF) that has a resolution compatible with the widely used Cα protein representation [D. Ugarte La Torre and S. Takada, J. Chem. Phys. 153, 205101 (2020)]. In this study, we extended it and developed a lipid-protein interaction model that allows the combination of the iSoLF and the Cα protein force field, AICG2+. The hydrophobic-hydrophilic interaction is modeled as a modified Lennard-Jones potential in which parameters were tuned partly to reproduce the experimental transfer free energy and partly based on the free energy profile normal to the membrane surface from previous all-atom MD simulations. Then, the obtained lipid-protein interaction is tested for the configuration and placement of transmembrane proteins, water-soluble proteins, and peripheral proteins, showing good agreement with prior knowledge. The interaction is generally applicable and is implemented in the publicly available software, CafeMol.

Stabilizing Leaflet Asymmetry in a Highly Coarse-Grained Lipid Membrane Model

Lipid bilayer asymmetry is an essential aspect of biomembranes that leads to properties such as spontaneous curvature and differential stress, which in turn influence the large-scale membrane geometry. This presents a dilemma in simulations: models at atomistic resolution require potentially prohibitive computational resources, while highly coarse-grained models struggle to maintain asymmetry due to accelerated flip-flop dynamics. To address this issue, we present here a version of the coarse-grained Cooke lipid model, amended to stabilize asymmetric lipid membranes.

Coarse-grained implicit solvent lipid force field with a compatible resolution to the Cα protein representation.

Biological membranes have been prominent targets for coarse-grained (CG) molecular dynamics simulations. While minimal CG lipid models with three beads per lipid and quantitative CG lipid models with >10 beads per lipid have been well studied, in between them, CG lipid models with a compatible resolution to residue-level CG protein models are much less developed. Here, we extended a previously developed three-bead lipid model into a five-bead model and parameterized it for two phospholipids, POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) and DPPC (1,2-dipalmitoyl-sn-glycero-3-phosphatidylcholine). The developed model, iSoLF, reproduced the area per lipid, hydrophobic thickness, and phase behaviors of the target phospholipid bilayer membranes at the physiological temperature. The model POPC and DPPC membranes were in liquid and gel phases, respectively, in accordance with experiments. We further examined the spontaneous formation of a membrane bilayer, the temperature dependence of physical properties, the vesicle dynamics, and the POPC/DPPC two-component membrane dynamics of the CG lipid model, showing some promise. Once combined with standard Cα protein models, the iSoLF model will be a powerful tool to simulate large biological membrane systems made of lipids and proteins.

Stabilizing Leaflet Asymmetry under Differential Stress in a Highly Coarse-Grained Lipid Membrane Model.

We present a version of the coarse-grained Cooke lipid model, modified to simulate asymmetric lipid membranes. It is inspired by a method employed by Wang et al. [ Commun. Comput. Phys. 2013, 13, 1093-1106] for artificially penalizing lipid flip-flop but copes more robustly with differential stress, at the cost of one additional bead per lipid and the concomitant increase in computational overhead. Bilayer asymmetry ultimately breaks down beyond a system size dependent critical differential stress, which can be predicted from a simple analytical model. We remeasure many important material parameters for the new model and find it to be consistent with typical fluid lipid membranes. Maintaining a stable stress asymmetry has many applications, and we give two examples: (i) connecting monolayer stress to lipid number asymmetry in order to directly measure the monolayer area modulus and (ii) finding its strain-dependent higher-order correction by monitoring the equilibrium bilayer area.

Development of a coarse-grained lipid model, LIME 2.0, for DSPE using multistate iterative Boltzmann inversion and discontinuous molecular dynamics simulations

We suggest an improved version of the intermediate resolution implicit solvent model for lipids, LIME, that was previously developed for use with discontinuous molecular dynamics (DMD) simulations. LIME gets its geometrical and energy parameters between bonded and nonbonded pairs of coarse-grained (CG) sites from atomistic simulations. The improved model, LIME 2.0, uses multiple square wells rather than the single square well used in original LIME to obtain intermolecular interactions that more faithfully mimic those from atomistic simulations. The multi-state iterative Boltzmann inversion (MS-IBI) scheme is used to determine the interaction parameters. This means that a single set of interaction parameters between coarse-grained sites can be used to represent the lipid bilayers at different temperatures. The physical properties of CG DSPE lipid bilayer are calculated using CG simulations and compared to atomistic simulations results to verify the improved model. The phase transition temperature of the lipid bilayer is measured accurately and the lipid translocation phenomenon, “flip-flop” is observed through CG simulation. These results suggest that CG parameterization using multiple square-wells and the MS-IBI scheme is well suited to the study of lipid bilayers across a range of temperatures with DMD simulations.

Transition States of Passive Lipid Transport are Characterized by Hydrophobic Contacts

Cell homeostasis requires the maintenance of heterogeneous membrane compositions through both vesicular and non-vesicular transport mechanisms. Yet, our biophysical description of non-vesicular lipid transport remains incomplete. To help fill these gaps, we aimed to identify an accurate reaction coordinate for passive lipid exchange between membranes. Towards this goal, we investigated the elementary steps of lipid exchange, lipid desorption and insertion into a membrane, using molecular simulation. From over 1,000 lipid insertion trajectories of all-atom and coarse-grained lipid models, we discovered a free energy barrier for lipid insertion and identified multiple pathways characterized by splayed lipid intermediates. This barrier appears hidden when only the lipid's displacement normal to the bilayer, which has traditionally been used to describe lipid exchange, is monitored. In contrast, an accurate reaction coordinate measures the breakage and formation of hydrophobic lipid-membrane contacts, which give rise to a barrier for lipid insertion. At the transition state, hydrophobic contacts are just as likely to form as they are to break. Consistent with this fact, membrane distortions and solvent fluctuations, which can both enable and prevent hydrophobic contact formation, are observed in the transition state ensemble. Overall, our results demonstrate that the formation and breakage of hydrophobic contacts is rate limiting for passive lipid exchange and provide a foundation to understand the catalytic function of lipid transfer proteins.

Extending a Highly Coarse-Grained Lipid Model to Asymmetric Membranes for MD Simulations

Lipid leaflet asymmetry is an essential aspect of biomembranes that leads to numerous novel properties not present in ordinary symmetric membranes. To investigate these systems through molecular dynamics simulation over large length and time scales usually necessitates coarse-grained models to reduce computational cost. However, highly coarse-grained models, such as the Cooke model, exhibit artificially accelerated dynamics, in particular lipid flip-flop. This may cause a membrane's asymmetry to decay before statistically significant measurements of its consequences can be made. We investigate a remedy for this problem proposed by Wang et al. (Commun. Comput. Phys. 13), which energetically discourages lipids from entering the opposite leaflet. We measure the efficiency of this modification in maintaining asymmetry in normal unstressed membranes as well as bilayers with varying degrees of stress imbalance between leaflets. We propose further changes to improve the remedy under extenuating circumstances, and we report changes to basic membrane characteristics caused by the modification as well as possible ways by which one can compensate for them.

Lipid flip-flop vs. lateral diffusion in the relaxation of hemifusion diaphragms

Molecular dynamics simulations of a solvent-free coarse-grained lipid model are used to characterize the mechanisms by which lipid-bilayer hemifusion diaphragm (HD) intermediates relax, across a range of global compositions of negative intrinsic curvature (NIC) lipids and neutral-curvature lipids. At low concentrations of NIC lipids, rapid fission produces a double bilayer end state through a lateral diffusion-based mechanism enabled by spontaneous rim-pore defects. At moderately higher NIC lipid concentrations, rim pores are absent and stable leaflet three-junctions persist, revealing an HD relaxation mechanism entirely reliant on lipid flip-flop, and end states that are either stable fusion pores or stable HD's. These fusogenic systems exhibit dynamics highly dependent on NIC lipid concentration via an underlying sensitivity of flip-flop rates for neutral lipids on NIC lipid concentration. This work illustrates that HD dynamics may be altered through regulation of lipid composition in the immediate three-junction region. This work further highlights the potential role of flippases in biological fusion and the importance of lipid composition on fusion dynamics.

Rate of hemifusion diaphragm dissipation and ability to form three-junction bound HD determined by lipid composition.

Though the hemifusion diaphragm (HD) is widely accepted as an intermediate in bilayer membrane fusion, lipid contributions toward HD stability and dynamics are still not fully understood. In this paper, we study large, binary, protein-free HD systems at varying compositions of negative intrinsic curvature (NIC) lipids using molecular dynamics (MD) simulations of a solvent-free coarse-grained lipid model. Under MD, initially created HDs are found to relax to three major end states depending on the composition and lipid intrinsic curvature. Low compositions of NIC lipids or weak intrinsic curvature result in double-bilayer end states, and moderate compositions of moderate to strong NIC lipids result in metastable fusion pores. Importantly, high compositions of moderate NIC lipids result in a metastable HD that persists beyond μs time scales. NIC lipids stabilize the HD by filling the junction core around the HD. Sorting of NIC lipids toward the three-junction region occurs in fused-endpoint systems, but no significant sorting was seen in systems that end in a double bilayer indicating that high line tension at the triple junction drives HD dissipation faster than sorting can enrich that junction enough to lower that line tension. The appearance of three end states dependent on the NIC lipid composition highlights the necessity of NIC lipids for non-leaky fusion.

Exploiting Lipid Permutation Symmetry to Compute Membrane Remodeling Free Energies.

A complete physical description of membrane remodeling processes, such as fusion or fission, requires knowledge of the underlying free energy landscapes, particularly in barrier regions involving collective shape changes, topological transitions, and high curvature, where Canham-Helfrich (CH) continuum descriptions may fail. To calculate these free energies using atomistic simulations, one must address not only the sampling problem due to high free energy barriers, but also an orthogonal sampling problem of combinatorial complexity stemming from the permutation symmetry of identical lipids. Here, we solve the combinatorial problem with a permutation reduction scheme to map a structural ensemble into a compact, nondegenerate subregion of configuration space, thereby permitting straightforward free energy calculations via umbrella sampling. We applied this approach, using a coarse-grained lipid model, to test the CH description of bending and found sharp increases in the bending modulus for curvature radii below 10nm. These deviations suggest that an anharmonic bending term may be required for CH models to give quantitative energetics of highly curved states.

Determining the Pivotal Plane of Fluid Lipid Membranes in Simulations

Each leaflet of a curved lipid membrane contains a surface at which the area strain vanishes, the so-called pivotal plane. Its distance z0 from the bilayer's mid-plane arises in numerous membrane structure related contexts, for instance the connection between monolayer and bilayer moduli, stress-profile moments, or area-difference elasticity theories. Experimentally, the pivotal plane is known to lie close to the glycerol backbone of a lipid, but only very few simulations have tried to locate it. Here we propose two precise methods for determining z0, both of which rely on monitoring the lipid imbalance across a curved bilayer. The first method considers the ratio of lipid number between the two leaflets of cylindrical or spherical vesicles and requires lipid flip-flop for equilibration. The second method looks at the leaflet difference across local sections cut out from a buckled membrane and is free of this limitation. Our simulations rely on two different coarse-grained lipid models. First, the generic three-bead solvent-free Cooke model, which is amenable to both methods and gives results for z0 that agree at the percent level. And second, a ten-bead representation of dimyristoylphosphocholine (DMPC) with the explicit solvent MARTINI model, which can only be analyzed with the buckling method. For the latter, the obtained value z0 = 0.850(11) nm lies about 0.4 nm inwards of the glycerol backbone, essentially in the middle of the leaflet, and is hence unexpectedly small. We attribute this to limitations of the coarse-grained description, suggesting that the location of the pivotal plane might be a novel indicator for how well lipid models capture the microscopic origins of curvature elasticity, especially the relative contributions of head and tail regions to the overall curvature modulus.

Determining the pivotal plane of fluid lipid membranes in simulations.

Each leaflet of a curved lipid membrane contains a surface at which the area strain vanishes, the so-called pivotal plane. Its distance z0 from the bilayer's midplane arises in numerous contexts, for instance the connection between monolayer and bilayer moduli, stress-profile moments, or area-difference elasticity theories. Here, we propose two precise methods for determining the location of the pivotal plane in computer simulations, both of which rely on monitoring the lipid imbalance across a curved bilayer. The first method considers the ratio of lipid number between the two leaflets of cylindrical or spherical vesicles; it hence requires lipid flip-flop for equilibration. The second method looks at the leaflet difference across local sections cut out from a buckled membrane; this observable equilibrates even in the absence of flip-flop. We apply our methods to two different coarse-grained lipid models, the generic three-bead solvent-free Cooke model and a ten-bead representation of dimyristoylphosphocholine with the explicit solvent MARTINI model. The Cooke model is amenable to both methods and gives results that agree at the percent level. Using it, we also show that the pivotal plane moves outward as lipid curvature becomes more positive. The MARTINI model can only be analyzed with the buckling method; the obtained value z0 = 0.850(11) nm lies about 0.4 nm inwards of the glycerol backbone and is hence unexpectedly small. We attribute this to limitations of the coarse-grained description, suggesting that the location of the pivotal plane might be a good indicator for how well lipid models capture the microscopic origins of curvature elasticity. Finally, we also show that the pivotal plane position itself moves as the membrane is bent. The leading correction is linear in curvature, dependent on the Poisson ratio, and can matter when analyzing experimental results obtained from highly curved inverse hexagonal phases.

Do Gel Phase Lipid Bilayers Behave Like Euler Elastica?

The elasticity of fluid phase lipid membranes can be characterized using the Helfrich Hamiltonian, which depends only on the square of the total curvature and the Gaussian curvature and the corresponding curvature moduli: the mean bending modulus and the Gaussian curvature modulus. Even at large curvatures approaching the inverse bilayer thickness, the effect of higher order terms of the elastic energy are minimal. Recently, a method has been developed to derive the bending modulus for fluid membranes from the stress-strain relationship of a buckled membrane [1,2]. The method also predicts the shape of the membrane, an Euler Elastica, and this serves as a check that the membrane indeed follows quadratic curvature elasticity. Using a coarse-grained lipid model, we analyze the shape of a buckled membrane in a gel phase and show that it does not behave like an Euler Elastica, even at low curvatures. The deviation from the theory suggests that higher order terms of the total curvature reduce the energy penalty for high curvatures. We present an extended version of the Helfrich Hamiltonian that captures this effect and show that it describes the shapes that we observe in simulations. We then calculate the bending modulus, as well as the modulus describing higher order corrections.[1] Noguchi H., “Anisotropic surface tension of buckled fluid membranes”, Phys. Rev. E 83, 061919 (2011).[2] Hu M., Diggins P., and Deserno M., “Determining the bending modulus of a lipid membrane by simulating buckling”, J. Chem. Phys. 138, 214110 (2013).

Mesoscale Simulation of Vesiculation of Lipid Droplets

An implicit solvent coarse-grained (CG) lipid model using three beads to reflect the basically molecular structure of two-tailed lipid is developed. In this model, the nonbonded interaction employs a variant MIE potential and the bonded interaction utilizes a Harmonic potential form. The CG force field parameters are achieved by matching the structural and mechanical properties of dipalmitoylphosphatidylcholine (DPPC) bilayers. The model successfully reproduces the formation of lipid bilayer from a random initial state and the spontaneous vesiculation of lipid bilayer from a disk-like structure. After that, the model is used to systematically study the vesiculation processes of spherical and cylindrical lipid droplets. The results show that the present CG model can effectively simulate the formation and evolution of mesoscale complex vesicles.

Solvent-Free, Highly Coarse-Grained Models for Charged Lipid Systems.

We present a methodology to develop coarse-grained lipid models such that electrostatic interactions between the coarse-grained sites can be derived accurately from an all-atom molecular dynamics trajectory and expressed as an effective pairwise electrostatic potential with appropriate screening functions. The reference nonbonded forces from the all-atom trajectory are decomposed into separate electrostatic and van der Waals (vdW) forces, based on the multiscale coarse-graining method. The coarse-grained electrostatic potential is assumed to be a general function of unknown variables and the final site–site interactions are obtained variationally, where the residual of the electrostatic forces from the assumed field is minimized. The resulting electrostatic interactions are fitted to screened electrostatics functions, with a special treatment for distance-dependent dielectrics and screened dipole–dipole interactions. The vdW interactions are derived separately. The resulting charged hybrid coarse-graining method is applied to various solvent-free three-site models of anionic lipid systems.

Buckling Gel-Phase Membranes is a Way to Measure their Mean Bending Regidity

The elasticity of lipid membranes can be characterized by two curvature moduli: the mean bending modulus and the Gaussian curvature modulus. Due to the relevance of the mean bending modulus for countless biological processes, considerable effort has been devoted to determine it_both in experiment and in computational studies. The most common computational approach is to measure the power spectrum of shape undulation modes [1]. Unfortunately, this technique is challenging for gel-phase membranes, because their larger modulus renders their fluctuations correspondingly smaller. In contrast, methods that infer a membrane's rigidity from actively bending it yield a signal that becomes stronger as the membrane becomes stiffer. One recently proposed method derives the modulus from the stress-strain relation of a buckled membrane, and it has been shown to provide accurate results for fluid membranes [2,3]. Using a coarse-grained lipid model, we show that this buckling method can also calculate the mean bending modulus of a gel phase membrane. We discuss the efficient implementation of the technique, paying special attention to difficulties that can arise while simulating a membrane in the gel phase. The method also provides insights into the contribution of entropy to the bending modulus, and hence its the temperature dependence.[1] Goetz R., Gompper G., and Lipowsky R., “Mobility and elasticity of self-assembled membranes”, Phys. Rev. Lett. 82, 221-224 (1999).[2] Noguchi H., “Anisotropic surface tension of buckled fluid membranes”, Phys. Rev. E 83, 061919 (2011).[3] Hu M., Diggins P., and Deserno M., “Determining the bending modulus of a lipid membrane by simulating buckling”, J. Chem. Phys. 138, 214110 (2013).

Embedding Aβ42 in Heterogeneous Membranes Depends on Cholesterol Asymmetries

Using a coarse-grained lipid and peptide model, we show that the free energy stabilization of amyloid-β in heterogeneous lipid membranes is predicted to have a dependence on asymmetric distributions of cholesterol compositions across the membrane leaflets. We find that a highly asymmetric cholesterol distribution that is depleted on the exofacial leaflet but enhanced on the cytofacial leaflet of the model lipid membrane thermodynamically favors membrane retention of a fully embedded Aβ peptide. However, in the case of cholesterol redistribution that increases concentration of cholesterol on the exofacial layer, typical of aging or Alzheimer’s disease, the free energy favors peptide extrusion of the highly reactive N-terminus into the extracellular space that may be vulnerable to aggregation, oligomerization, or deleterious oxidative reactivity.