Simulations Of Lipid Bilayers Research Articles

Simulating asymmetric membranes using P21 periodic boundary conditions.

Molecular dynamics (MD) simulations of symmetric lipid bilayers are now well established, while those of asymmetric ones are considerably less developed. This disjunction arises in part because the surface tensions of leaflets in asymmetric bilayers can differ (unlike those of symmetric ones), and there is no simple way to determine them without assumptions. This chapter describes the use of P21 periodic boundary conditions (PBC), which allow lipids to switch leaflets, to generate asymmetric bilayers under the assumption of equal chemical potentials of lipids in opposing leaflets. A series of examples, ranging from bilayers with one lipid type to those with peptides and proteins, provides a guide for the use of P21 PBC. Critical properties of asymmetric membranes, such as spontaneous curvature, are highly sensitive to differences in the leaflet surface tensions (or differential stress), and equilibration with P21 PBC substantially reduces differential stress of asymmetric bilayers assembled with surface area-based methods. Limitations of the method are discussed. Technically, the nonstandard unit cell is difficult to parallelize and to incorporate restraints. Inherently, the assumption of equal chemical potentials, and therefore the method itself, is not applicable to all target systems. Despite these limitations, it is argued that P21 simulations should be considered when designing equilibration protocols for MD studies of most asymmetric membranes.

Dynamics behavior of PE and PET oligomers in lipid bilayer simulations

In recent years many investigators have been concerned about the toxicity and potential health hazards of micro- and nanoplastics. However, we are still lacking a good understanding of the methods of their transport into the human body and subsequently within cells. This is especially true at the lower nanometer scale; these particles are potentially more dangerous than their micrometer counterparts due to their easier permeation into cells. In this study we used both unbiased molecular dynamics simulations and steered umbrella sampling simulations to explore the interactions of polyethylene terephthalate (PET) and polyethylene (PE) oligomers in phospholipid bilayers. Our simulations revealed that the bilayers did not represent significant energy barriers to the small oligomers; not only did they readily enter the cell membrane but they also became concentrated into specific parts of the membrane. The larger PET tetramers exhibited a strong aggregation in water but were the least likely to permeate through or into the membranes. It is possible that PE monomers and tetramers can become concentrated into membranes while PET monomers are more likely to pass through or concentrate just inside the membrane surface. Passive transport of microplastics into cells is, however, likely limited to particles of a few nanometers in diameter.

Understanding the Mechanical Properties of Ultradeformable Liposomes Using Molecular Dynamics Simulations.

Improving drug delivery efficiency to solid tumor sites is a central challenge in anticancer therapeutic research. Our previous experimental study (Guo et al., Nat. Commun. 2018, 9, 130) showed that soft, elastic liposomes had increased uptake and accumulation in cancer cells and tumors in vitro and in vivo respectively, relative to rigid particles. As a first step toward understanding how liposomes' molecular structure and composition modulates their elasticity, we performed all-atom and coarse-grained classical molecular dynamics (MD) simulations of lipid bilayers formed by mixing a long-tailed unsaturated phospholipid with a short-tailed saturated lipid with the same headgroup. The former types of phospholipids considered were 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) and 1,2-dipalmitoleoyl-sn-glycero-3-phosphocholine (termed here DPMPC). The shorter saturated lipids examined were 1,2-diheptanoyl-sn-glycero-3-phosphocholine (DHPC), 1,2-didecanoyl-sn-glycero-3-phosphocholine (DDPC), 1,2-dilauroyl-sn-glycero-3-phosphocholine (DLPC), and 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC). Several lipid concentrations and surface tensions were considered. Our results show that DOPC or DPMPC systems having 25-35 mol % of the shortest lipids DHPC or DDPC are the least rigid, having area compressibility moduli KA that are ∼10% smaller than the values observed in pure DOPC or DPMPC bilayers. These results agree with experimental measurements of the stretching modulus and lysis tension in liposomes with the same compositions. These mixed systems also have lower areas per lipid and form more uneven x-y interfaces with water, the tails of both primary and secondary lipids are more disordered, and the terminal methyl groups in the tails of the long lipid DOPC or DPMPC wriggle more in the vertical direction, compared to pure DOPC or DPMPC bilayers or their mixtures with the longer saturated lipid DLPC or DMPC. These observations confirm our hypothesis that adding increasing concentrations of the short unsaturated lipid DHPC or DDPC to DOPC or DPMPC bilayers alters lipid packing and thus makes the resulting liposomes more elastic and less rigid. No formation of lipid nanodomains was noted in our simulations, and no clear trends were observed in the lateral diffusivities of the lipids as the concentration, type of secondary lipid, and surface tension were varied.

Cellular excitability and ns-pulsed electric fields: Potential involvement of lipid oxidation in the action potential activation

Recent studies showed that nanosecond pulsed electric fields (nsPEFs) can activate voltage-gated ion channels (VGICs) and trigger action potentials (APs) in excitable cells. Under physiological conditions, VGICs’ activation takes place on time scales of the order 10–100 µs. These time scales are considerably longer than the applied pulse duration, thus activation of VGICs by nsPEFs remains puzzling and there is no clear consensus on the mechanisms involved. Here we propose that changes in local electrical properties of the cell membrane due to lipid oxidation might be implicated in AP activation. We first use MD simulations of model lipid bilayers with increasing concentration of primary and secondary lipid oxidation products and demonstrate that oxidation not only increases the bilayer conductance, but also the bilayer capacitance. Equipped with MD-based characterization of electrical properties of oxidized bilayers, we then resort to AP modelling at the cell level with Hodgkin-Huxley-type models. We confirm that a local change in membrane properties, particularly the increase in membrane conductance, due to formation of oxidized membrane lesions can be high enough to trigger an AP, even when no external stimulus is applied. However, excessive accumulation of oxidized lesions (or other conductive defects) can lead to altered cell excitability.

Extension of the iSoLF implicit-solvent coarse-grained model for multicomponent lipid bilayers.

iSoLF is a coarse-grained (CG) model for lipid molecules with the implicit-solvent approximation used in molecular dynamics (MD) simulations of biological membranes. Using the original iSoLF (iSoLFv1), MD simulations of lipid bilayers consisting of either POPC or DPPC and these bilayers, including membrane proteins, can be performed. Here, we improve the original model, explicitly treating the electrostatic interactions between different lipid molecules and adding CG particle types. As a result, the available lipid types increase to 30. To parameterize the potential functions of the new model, we performed all-atom MD simulations of each lipid at three different temperatures using the CHARMM36 force field and the modified TIP3P model. Then, we parameterized both the bonded and non-bonded interactions to fit the area per lipid and the membrane thickness of each lipid bilayer by using the multistate Boltzmann Inversion method. The final model reproduces the area per lipid and the membrane thickness of each lipid bilayer at the three temperatures. We also examined the applicability of the new model, iSoLFv2, to simulate the phase behaviors of mixtures of DOPC and DPPC at different concentrations. The simulation results with iSoLFv2 are consistent with those using Dry Martini and Martini 3, although iSoLFv2 requires much fewer computations. iSoLFv2 has been implemented in the GENESIS MD software and is publicly available.

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.

Surface viscosities of several lipids determined from equilibrium molecular dynamics simulations.

Lipid membrane viscosity is critical to biological function. Bacterial cells grown in different environments have been shown to alter their lipid composition in order to maintain a viscosity in a narrow range,1 and membrane viscosity has been shown to control the rate of cellular respiration.2 In order to understand the factors which determine membrane viscosity, 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 fit to 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 twenty times that of DPPC. Including the full range of the dispersion increases the membrane viscosity by up to 140%, alleviating some of the discrepancy between simulated and experimental viscosity measurements, although the simulated viscosities remain significantly lower than recent experimental values reported for a series of lipids.3

Efficient calculation of the free energy for protein partitioning using restraining potentials

An approach for the efficient simulation of phase-separated lipid bilayers, for use in the calculation of equilibrium free energies of partitioning between lipid domains, is proposed. The methodology exploits restraint potentials and rectangular aspect ratios that enforce lipid phase separation, allowing for the simulation of smaller systems that approximately reproduce bulk behavior. The utility of this approach is demonstrated through the calculation of potentials of mean force for the translation of a transmembrane protein between lipid domains. The impact of the imposed restraints on lipid tail ordering and lipid packing are explored, providing insight into how restraints can best be employed to compute accurate free-energy surfaces. This approach should be useful in the accurate calculation of equilibrium partition coefficients for transmembrane protein partitioning in heterogeneous membranes, providing insight into the thermodynamic driving forces that control this fundamental biophysical phenomenon.

Molecular substructure of the liquid-ordered phase formed by sphingomyelin and cholesterol: sphingomyelin clusters forming nano-subdomains are a characteristic feature.

As a model of lipid rafts, the liquid-ordered (Lo) phase formed by sphingomyelin (SM) and cholesterol (Cho) in bilayer membranes has long attracted the attention of biophysics researchers. New approaches and methodologies have led to a better understanding of the molecular basis of the Lo domain structure. This review summarizes studies on model membrane systems consisting of SM/unsaturated phospholipid/Cho implying that the Lo phase contains SM-based nanodomains (or nano-subdomains). Some of the Lo phase properties may be attributed to these nanodomains. Several studies suggest that the nanodomains contain clustered SM molecules packed densely to form gel-phase-like subdomains of single-digit nanometer size at physiological temperatures. Cho and unsaturated lipids located in the Lo phase are likely to be concentrated at the boundaries between the subdomains. These subdomains are not readily detected in the Lo phase formed by saturated phosphatidylcholine (PC) molecules, suggesting that they are strongly stabilized by homophilic interactions specific to SM, e.g., between SM amide groups. This model for the Lo phase is supported by experiments using dihydro-SM, which is thought to have stronger homophilic interactions than SM, as well as by studies using the enantiomer of SM having opposite stereochemistry to SM at the 2 and 3 positions and by some molecular dynamics (MD) simulations of lipid bilayers containing Lo-lipids. Nanosized gel subdomains seem to play an important role in controlling membrane organization and function in biological membranes.

Algorithm to catalogue topologies of dynamic lipid hydrogen-bond networks.

Lipid membrane interfaces host reactions essential for the functioning of cells. The hydrogen-bonding environment at the membrane interface is particularly important for binding of proteins, drug molecules, and ions. We present here the implementation and applications of a depth-first search algorithm that analyzes dynamic lipid interaction networks. Lipid hydrogen-bond networks sampled transiently during simulations of lipid bilayers are clustered according to main types of topologies that characterize three-dimensional arrangements of lipids connected to each other via short water bridges. We characterize the dynamics of hydrogen-bonded lipid clusters in simulations of model POPE and POPE:POPG membranes that are often used for bacterial membrane proteins, in a model of the Escherichia coli membrane with six different lipid types, and in POPS membranes. We find that all lipids sample dynamic hydrogen-bonded networks with linear, star, or circular arrangements of the lipid headgroups, and larger networks with combinations of these three types of topologies. Overall, linear lipid-water bridges tend to be short. Water-mediated lipid clusters in all membranes with PE lipids tend to be somewhat small, with about four lipids in all membranes studied here. POPS membranes allow circular arrangements of three POPS lipids to be sampled frequently, and complex arrangements of linear, star, and circular paths may also be sampled. These findings suggest a molecular picture of the membrane interface whereby lipid molecules transiently connect in clusters with somewhat small spatial extension.

WIND-PVPA: Water/Ion NMR Detected PVPA to assess lipid barrier integrity in vitro through quantification of passive water- and ion transport

Water/Ion NMR Detected – Phospholipid Vesicle Permeability Assay (WIND-PVPA), is presented as a novel, straightforward and automatable method to assess lipid barrier integrity in vitro. The apparent permeability constants of water- and ions across the PVPA barriers are determined in a one-pot experiment under the influence of membrane-active guest molecules. NMR spectroscopy is used to quantify the water directly (D2O) and the ions indirectly (complexed with EDTA) as a function of time. WIND-PVPA is demonstrated using four anti-microbial peptides, to show that membrane active molecules can be differentiated by their disruptive influence on the PVPA system. The results obtained are compared with explicit molecular dynamics simulations of lipid bilayers, AMPs, water and salt, where the motions of all individual water molecules relative to the lipid bilayer are monitored over the course of the simulations, allowing the calculation of theoretical apparent permeability constants of the corresponding single bilayer systems.Proof-of-principle is presented that WIND-PVPA can be used to evaluate the lipid barrier destabilizing effect of active guest molecules by measuring changes in passive water- and ion permeabilities upon exposure. The method is highly flexible in terms of barrier composition, choice of probes and membrane active compounds.

A Method for Detection of Water Permeation Events in Molecular Dynamics Simulations of Lipid Bilayers

A Method for Detection of Water Permeation Events in Molecular Dynamics Simulations of Lipid Bilayers

High-Dimensional Parameter Search Method to Determine Force Field Mixing Terms in Molecular Simulations.

Molecular dynamics (MD) force fields for lipids and ions are typically developed independently of one another. In simulations consisting of both lipids and ions, lipid-ion interaction energies are estimated using a predefined set of mixing rules for Lennard-Jones (LJ) interactions. This, however, does not guarantee their reliability. In fact, compared to the quantum mechanical reference data, Lorentz-Berthelot mixing rules substantially underestimate the binding energies of Na+ ions with small-molecule analogues of lipid headgroups, yielding errors on the order of 80 and 130 kJ/mol, respectively, for methyl acetate and diethyl phosphate. Previously, errors associated with mixing force fields have been reduced using approaches such as "NB-fix" in which LJ interactions are computed using explicit cross terms rather than those from mixing rules. Building on this idea, we derive explicit lipid-ion cross terms that also may implicitly include many-body cooperativity effects. Additionally, to account for the interdependency between cross terms, we optimize all cross terms simultaneously by performing high-dimensional searches using our ParOpt software. The cross terms we obtain reduce the errors due to mixing rules to below 10 kJ/mol. MD simulation of the lipid bilayer conducted using these optimized cross terms resolves the structural discrepancies between our previous simulations and small-angle X-ray and neutron scattering experiments. These results demonstrate that simulations of lipid bilayers with ions that are accurate up to structural data from scattering experiments can be performed without explicit polarization terms. However, it is worth noting that such NB-fix cross terms are not based on any physical principle; a polarizable lipid model would be more realistic and is still desired. Our approach is generic and can be applied to improve the accuracies of simulations employing mixed force fields.

The molecular determinants of membrane viscosity

Lipid membrane viscosity is critical to biological function. Bacterial cells grown in different environments have been shown to alter their lipid composition in order to maintain a specific viscosity, and membrane viscosity has been linked to the rate of cellular respiration. In order to understand the factors which 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 fit to a stretched exponential by a maximum-likelihood Markov chain Monte Carlo method.

The effect of different force fields on the ion-membrane interactions

Molecular dynamics (MD) simulations of lipid membranes have been widely used to better understand phenomena at the atomic and molecular level that cannot be observed with conventional experimental methods. This technique is very useful to understand the underlying mechanisms of the most common biological ions, such as K+, Na+, Cl- or Ca2+ when they interact with the lipid bilayer or when they go through the induced pores during electroporation. Lipid bilayer simulations showed that K+ and Cl- interact very little with the lipids, compared to Na+ and Ca2+.

Ripple-like instability in the simulated gel phase of finite size phosphocholine bilayers.

Atomistic molecular dynamics simulations have reached a degree of maturity that makes it possible to investigate the lipid polymorphism of model bilayers over a wide range of temperatures. However if both the fluid Lα and tilted gel [Formula: see text] states are routinely obtained, the [Formula: see text] ripple phase of phosphatidylcholine lipid bilayers is still unsatifactorily described. Performing simulations of lipid bilayers made of different numbers of DPPC (1,2-dipalmitoylphosphatidylcholine) molecules ranging from 32 to 512, we demonstrate that the tilted gel phase [Formula: see text] expected below the pretransition cannot be obtained for large systems (equal or larger than 94 DPPC molecules) through common simulations settings or temperature treatments. Large systems are instead found in a disordered gel phase which display configurations, topography and energies reminiscent from the ripple phase [Formula: see text] observed between the pretransition and the main melting transition. We show how the state of the bilayers below the melting transition can be controlled and depends on thermal history and conditions of preparations. A mechanism for the observed topographic instability is suggested.

Nonconverged Constraints Cause Artificial Temperature Gradients in Lipid Bilayer Simulations.

Molecular dynamics (MD) simulations have become an indispensable tool to investigate phase separation in model membrane systems. In particular, simulations based on coarse-grained (CG) models have found widespread use due to their increased computational efficiency, allowing for simulations of multicomponent lipid bilayers undergoing phase separation into liquid-ordered and liquid-disordered domains. Here, we show that a significant temperature difference between molecule types can artificially arise in CG MD membrane simulations with the standard Martini simulation parameters in GROMACS. In particular, the linear constraint solver (LINCS) algorithm does not converge with its default settings, resulting in serious temperature differences between molecules in a time step-dependent manner. We demonstrate that the underlying reason for this behavior is the presence of highly constrained moieties, such as cholesterol. Their presence can critically impact numerous structural and dynamic membrane properties obtained from such simulations. Furthermore, any preference of these molecules toward a certain membrane phase can lead to spatial temperature gradients, which can amplify the degree of phase separation or even induce it in compositions that would otherwise mix well. We systematically investigated the effect of the integration time step and LINCS settings on membrane properties. Our data show that for cholesterol-containing membranes, a time step of 20 fs should be combined with at least lincs_iter = 2 and lincs_order = 12, while using a time step of 30 fs requires at least lincs_iter = 3 and lincs_order = 12 to bring the temperature differences to a level where they do not perturb central membrane properties. Moreover, we show that in cases where stricter LINCS settings are computationally too demanding, coupling the lipids in multiple groups to the temperature bath offers a practical workaround to the problem, although the validity of this approach should be further verified. Finally, we show that similar temperature gradients can also emerge in atomistic simulations using the CHARMM force field in combination with settings that allow for a 5 fs integration step.

Impact of Increased Membrane Realism on Conformational Sampling of Proteins.

The realism and accuracy of lipid bilayer simulations through molecular dynamics (MD) are heavily dependent on the lipid composition. While the field is pushing toward implementing more heterogeneous and realistic membrane compositions, a lack of high-resolution lipidomic data prevents some membrane protein systems from being modeled with the highest level of realism. Given the additional diversity of real-world cellular membranes and protein-lipid interactions, it is still not fully understood how altering membrane complexity affects modeled membrane protein functions or if it matters over long-timescale simulations. This is especially true for organisms whose membrane environments have little to no computational study, such as the plant plasma membrane. Tackling these issues in tandem, a generalized, realistic, and asymmetric plant plasma membrane with more than 10 different lipid species is constructed herein. Classical MD simulations of pure membrane constructs were performed to evaluate how altering the compositional complexity of the membrane impacted the plant membrane properties. The apo form of a plant sugar transporter, OsSWEET2b, was inserted into membrane models where lipid diversity was calculated in either a size-dependent or size-independent manner. An adaptive sampling simulation regime validated by Markov-state models was performed to capture the gating dynamics of OsSWEET2b in each of these membrane constructs. In comparison to previous OsSWEET2b simulations performed in a pure POPC bilayer, we confirm that simulations performed within a native-like membrane composition alter the stabilization of apo OsSWEET2b conformational states by ∼1 kcal/mol. The free-energy barriers of intermediate conformational states decrease when realistic membrane complexity is simplified, albeit roughly within sampling error, suggesting that protein-specific responses to membranes differ due to altered packing caused by compositional fluctuations. This work serves as a case study where a more realistic bilayer composition makes unbiased conformational sampling easier to achieve than with simplified bilayers.

Simulations of Diabetic and Non-Diabetic Peripheral Nerve Myelin Lipid Bilayers

The multilayered myelin sheath is a critical component of both central and peripheral nervous systems, forming a protective barrier against axonal damage and facilitating the movement of nervous impulses. It is primarily composed of cholesterol (CHL1), phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylserine (PS), phosphatidylinositol (PI), sphingomyelin (SM), and galactosylceramide (GalCer) lipids. For rat sciatic nerve myelin (part of the peripheral nervous system, PNS), it has been found that cholesterol and unsaturated fatty acid contents are significantly lower in diabetic than in non-diabetic conditions. In this study, lipid compositions from experimental data are used to create four model rat sciatic nerve myelin lipid bilayers: PI-containing (non-diabetic and diabetic) and PS-containing (non-diabetic and diabetic), which were then simulated using the all-atom CHARMM36 force field. Simulation results of diabetic membranes indicate less rigid, more laterally expansive, and thinner bilayers as well as potentially reduced interactions between GalCer on opposing myelin leaflets, supporting a direct role of the cholesterol content decrease in instigating myelin deterioration and diabetic peripheral neuropathy. Compared to PI-lipids, PS-lipids were found to cause higher inter-lipid spacing and decreased order within membranes as a result of their smaller headgroup size and higher inter-lipid hydrogen bonding potential, which allow them to more frequently reside deeper in the membrane plane and produce pushing effects on other lipids. GalCer deuterium order parameters and non-diabetic headgroup-to-headgroup bilayer thicknesses were compared to experimental data, exhibiting close alignment, supporting the future usage of these models to study the PNS myelin sheath.

The Extended Eighth-Shell method for periodic boundary conditions with rotational symmetry.

The Eighth-Shell method for parallelization of molecular dynamics simulations has previously been shown to be the most optimal for efficiency at large process counts. However, in its current formulation only the P1 space group is supported for periodic boundary conditions (PBC) and thus reflection and/or rotational crystal symmetries are not supported. In this work, we outline the development and implementation of the Extended Eighth-Shell (EES) method that allows rotational symmetry by using an extended import region compared to the ES method. It simulates only the asymmetric unit and communicates coordinates and forces with images that correspond to P21 PBC. The P21 PBC has application in lipid bilayer simulations as it can be used to allow lipids to switch leaftlets, thus rapidly balancing the chemical potential difference between the two layers. Our results show that the EES method scales efficiently over large number of processes and can be used for simulations with P21 symmetry in an orthorhombic crystal.