Micelle cluster distributions from molecular dynamics simulations of a solvent-free coarse-grained model of sodium octyl sulfate (SOS) were analyzed using an improved method to extract equilibrium association constants from small-system simulations containing one or two micelle clusters at equilibrium with free surfactants and counterions. The statistical–thermodynamic and mathematical foundations of this partition-enabled analysis of cluster histograms (PEACH) approach are presented. A dramatic reduction in computational time for analysis was achieved through a strategy similar to the selector variable method to circumvent the need for exhaustive enumeration of the possible partitions of surfactants and counterions into clusters. Using statistics from a set of small-system (up to 60 SOS molecules) simulations as input, equilibrium association constants for micelle clusters were obtained as a function of both number of surfactants and number of associated counterions through a global fitting procedure. The resulting free energies were able to accurately predict micelle size and charge distributions in a large (560 molecule) system. The evolution of micelle size and charge with SOS concentration as predicted by the PEACH-derived free energies and by a phenomenological four-parameter model fit, along with the sensitivity of these predictions to variations in cluster definitions, are analyzed and discussed.

Simulations of small unilamellar lipid bilayer vesicles have been performed to model their response to an instantaneous rise in temperature, starting from an initial low-temperature structure, to temperatures near or above the main chain transition temperature. The MARTINI coarse-grained force-field was used to construct slabs of gel-phase DPPC bilayers, which were assembled into truncated icosahedral structures containing 13165 or 31021 lipids. Equilibration at 280 K produced structures with several (5–8) domains, characterized by facets of lipids packed in the gel phase connected by disordered ridges. Instantaneous heating to final temperatures ranging from 290 K to 310 K led to partial or total melting over 500 ns trajectories, accompanied by changes in vesicle shape and the sizes and arrangements of remaining gel-phase domains. At temperatures that produced partial melting, the gel-phase lipid content of the vesicles followed an exponential decay, similar in form and timescale to the sub-microsecond phase of melting kinetics observed in recent ultrafast IR temperature-jump experiments. The changing rate of melting appears to be the outcome of a number of competing contributions, but changes in curvature stress arising from the expansion of the bilayer area upon melting are a major factor. The simulations give a more detailed picture of the changes that occur in frozen vesicles following a temperature jump, which will be of use for the interpretation of temperature-jump experiments on vesicles.

The umbrella sampling method has been used to evaluate the free energy profile for a large permeant moving through a lipid bilayer, represented using a coarse-grained simulation model, at and below its gel–fluid transition temperature. At the lipid transition temperature, determined to be 302 K for the MARTINI 2.0 model of DPPC, the permeation barrier for passage through an enclosed fluid domain embedded in a patch of gel was significantly lower than that for passage through a fluid stripe domain. In contrast, permeation through a fluid domain in a stripe geometry produced a free energy profile nearly identical to that of a gel-free fluid bilayer. In both cases, insertion of the permeant into a fluid domain coexisting with the gel phase led to a shift in phase composition, as lipids transitioned from fluid to gel to accommodate the area occupied by the permeant. In the case of the enclosed fluid domain, this transition produced a decrease in the length of the fluid–gel interface as the approximately circular fluid domain shrank. The observed decrease in the apparent permeation barrier, combined with an approximation for the change in interfacial length, enabled estimation of the interfacial line tension to be between 10 and 13 pN for this model. The permeation barrier was shown to drop even further in simulations performed at temperatures below the transition temperature. The results suggest a mechanism to explain the experimentally observed anomalous peak in the temperature-dependent permeability of lipid bilayers near their transition temperatures. The contribution of this mechanism toward the permeability of a gel phase containing a thermal distribution of fluid-phase domains is estimated using a simple statistical thermodynamic model.

When a range of lipid bilayers are melted to the disordered fluid phase from the (much less permeable) ordered gel phase, their permeability to a variety of permeants shows a peak at the transition temperature and drops off with increasing temperature, rather than just rising as melting proceeds. To explore this anomalous behavior, a simulated coarse-grained lipid membrane model that exhibits a phase transition upon expansion or compression was studied to determine how the permeation rate of a simple particle depends on the phase composition in the two-phase region and on particle size. The permeation rate and each phase’s area fraction and area density could be directly calculated, along with the probability that the permeant would cross in either phase or in interfacial regions. For large permeants and system sizes, conditions could be found where permeability increases upon compression of the bilayer. Permeation was negligible in the gel phase and, in contrast to the predictions of the “leaky interface” hypothesis, was not enriched in interfacial regions. The anomalous effect could instead be attributed to an increase in the area per lipid of fluid-phase domains. This result motivated a model for the decrease in effective permeability barrier through fluid-phase domains arising from a decrease in the length of the gel/fluid interface at the midpoint of a permeation event.

The structural properties and thermal stability of dipalmitoylphosphatidylethanolamine (DPPE) in the ordered gel phase have been studied by molecular dynamics simulation using two force fields: the Berger united-atom model and the CHARMM C36 atomistic model. As is widely known, structural features are sensitive to the initial preparation of the gel phase structure, as some degrees of freedom are slow to equilibrate on the simulation time scale of hundreds of nanoseconds. In particular, we find that the degree of alignment of the lipids’ glycerol backbones, which join the two hydrocarbon tails of each molecule, strongly affects the tilt angle of the tails in the resulting structures. Disorder in the backbone correlates with lower tilt angles: bilayer configurations initiated with aligned backbones produced tilt angles near 21° and 29° for the Berger and C36 force fields, respectively, while structures initiated with randomized backbone orientations showed average tilt angles of 7° and 18°, in closer agreement with the untilted structure observed experimentally. The transition temperature for the Berger force field gel bilayer has been determined by monitoring changes in width of gel phase stripe domains as a function of temperature and is 12 ± 5 K lower than the experimental value.

An approach is given to analyze aggregate size distributions obtained from simulations of a fixed number N of monomers undergoing reversible self-assembly. Equilibrium distributions are derived from size-dependent equilibrium association constants by appropriately weighted sums over all partitions of N monomers into aggregates. Conversely, equilibrium association constants can be obtained from an iterative fit to a finite-N equilibrium distribution. Model data for a micelle-forming system are used to show how results from simulations containing few micelles can yield infinite-N limiting distributions. A strategy is also suggested to exploit small-N effects on aggregate size distributions to enhance sampling of critical clusters in determination of nucleation free energy functions.

The long-time reorganization of peptide networks formed by mixing “junction” and “linker” components, where the junctions contain multiple high-affinity, high-specificity binding sites for divalent “linkers”, is investigated using a mesoscale simulation model. An interesting feature of this type of system is that the concentration of defects can be controlled experimentally by varying the ratio of junction and linker concentrations in the system. In this study we use a simple simulation model to evaluate how the migration of defects controls the rate of relaxation of shear stress, and how this rate is related to junction multiplicity, defect concentration, and linker stiffness. The mean stress relaxation per defect migration event obtained through simulation was two to three times greater than assumed by standard simple theories. The results are used to develop a phenomenological theory to predict the dependence of shear viscosity on the equilibrium constant and rate constants associated with linker-junction binding, stoichiometric mismatch, and network topology. The shear stress relaxation time and viscosity are predicted to fall away sharply when junction and linker concentrations are mismatched by even 1%. Furthermore, the time-dependence of viscosity during gel “aging” is modeled as an approach to equilibrium through diffusion-limited recombination of complementary defects.

Mixed MD/MC simulation at fixed difference in chemical potential (Δμ) between two lipid types provides a computational indicator of the relative affinities of the two lipids for different environments. Applying this technique to ternary DPPC/DOPC/cholesterol bilayers yields a DPPC/DOPC ratio that increases with increasing cholesterol content at fixed Δμ, consistent with the known enrichment of DPPC and cholesterol-rich in liquid-ordered phase domains in the fluid−fluid coexistence region of the ternary phase diagram. Comparison of the cholesterol-dependence of PC compositions at constant Δμ with experimentally measured coexistence tie line end point compositions affords a direct test of the faithfulness of the atomistic model to experimental phase behavior. DPPC/DOPC ratios show little or no dependence on cholesterol content at or below 16% cholesterol in the DOPC-rich region of the composition diagram, indicating cooperativity in the favorable interaction between DPPC and cholesterol. The relative affinity of DPPC and DOPC for high cholesterol bilayer environments in simulations is explicitly shown to depend on the degree of cholesterol alignment with the bilayer normal, suggesting that a source of the cooperativity is the composition dependence of cholesterol tilt angle distributions.

The local lipid composition near a transmembrane helical peptide in mixed-lipid bilayers has been studied using a mixed molecular dynamics (MD) and configuration-bias Monte Carlo method that allows the lateral distribution of lipids to equilibrate much more quickly than is possible by diffusive mixing alone. Gramicidin-A peptide was embedded in bilayer mixtures of DMPC with either DDPC (shorter by four carbons per tail) or DSPC (longer by four carbons per tail) at 330 K to investigate the possibility of lipid sorting by tail length. Conventional MD simulations showed local thickening of the bilayer near the peptide in pure DDPC and local thinning of the bilayer in pure DSPC, with comparatively little perturbation to the thickness of pure DMPC bilayers, suggesting that DMPC has the best matched tail length to the peptide of these three. In 1:1 DMPC:DDPC mixtures, the DMPC lipid was weakly enriched (by about 5%) near the peptide, while in DSPC:DMPC mixtures, no consistent trend was observed. The results underscore the weakness of the coupling between membrane deformation and local composition fluctuations in the absence of spontaneous phase separation.

Many lipid bilayers undergo a reversible order−disorder transition between the gel and liquid crystalline (LC) phases at a main phase transition temperature Tm that is an important characteristic property of the lipid. Although Tm should serve as a useful standard for validation and calibration of simulation models of lipid bilayers, its evaluation within simulations is difficult due to the slow kinetics of the gel−LC transition, especially near Tm. A stripe growth strategy for calculating Tm, which aims to bypass the slowest steps in this transition, has been applied to dipalmitoylphosphatidylcholine and distearoylphosphatidylcholine bilayers represented with a commonly used united-atom force field. The strategy consists of initial preparation of a bilayer containing gel and LC domains arranged as parallel stripes, observation of the direction and rate of domain growth over a range of temperatures, and fitting rates to an Arrhenius-like equation for their temperature dependence that crosses zero at Tm. Calculated Tm’s for both lipids are 5−6 degrees lower than their experimental values, in much closer agreement with experiment than suggested by recent simulations that simulate heating and cooling of bilayer patches. The stripe growth method also yields rates of phase front propagation that are in order-of-magnitude agreement with experimental estimates, as well as insight into glycerol backbone disordering at the LC−gel interface.

The slow rate of diffusive mixing poses a challenge for molecular dynamics (MD) simulation studies of mixed-lipid bilayers. A mixed Monte Carlo−molecular dynamics (MC−MD) approach, which uses mutation moves to swap lipid types throughout the system within the semi-grand canonical ensemble, is here applied to a comparison of binary mixtures in the gel and liquid crystalline phases. The two lipid components modeled, distearoylphosphatidylcholine (DSPC) and dimyristoylphosphatidylcholine (DMPC), differ by four carbons in the lengths of their acyl tails and are investigated here at full hydration at a temperature (313 K) between their transition temperatures, where coexistence between a DSPC-rich gel phase and a DMPC-rich liquid crystalline phase is expected. An analysis of DSPC-DMPC mixtures in the gel phase indicates strong deviation from ideality in the thermodynamics of mixing, accompanied by a tendency of the shorter-tailed component DMPC to associate laterally and for DMPC headgroups to be displaced toward the bilayer midplane. The liquid crystal phase mixtures, in contrast, show more mild deviation from thermodynamically ideal mixing with no apparent tendency for similar lipids to cluster laterally and no difference in headgroup normal distribution profiles.

To understand the effects of lipid composition on membrane protein function in a mixture as complex as a biomembrane, one must know whether the lipid composition local to the protein differs from the mean lipid composition. In this study, we simulated the transmembrane domain of a β-barrel protein, OmpA, in mixtures of lipids of different tail lengths under conditions of negative hydrophobic mismatch, i.e., local bilayer thinning. We modeled the influence of OmpA on the local lipid composition both at a coarse-grained (CG) resolution using conventional molecular dynamics, and at an atomistic resolution within the semi-grand canonical ensemble using mutation moves to rapidly approach an equilibrium lateral distribution of lipids. Moderate enrichment of the shorter tail component (either DDPC in DDPC/DMPC mixtures or DMPC in DMPC/DSPC mixtures) extending 2–3 nm away from the protein surface was observed with both the atomistic and CG models. The similarity in trends suggests that the more computationally economical CG models capture the essential features of lipid sorting induced by hydrophobic mismatch.

Mixtures of long- and short-tail phosphatidylcholine lipids are known to self-assemble into a variety of aggregates combining flat bilayerlike and curved micellelike features, commonly called bicelles. Atomistic simulations of bilayer ribbons and perforated bilayers containing dimyristoylphosphatidylcholine (DMPC, di-C14 tails) and dihexanoylphosphatidylcholine (DHPC, di-C6 tails) have been carried out to investigate the partitioning of these components between flat and curved microenvironments and the stabilization of the bilayer edge by DHPC. To approach equilibrium partitioning of lipids on an achievable simulation timescale, configuration-bias Monte Carlo mutation moves were used to allow individual lipids to change tail length within a semigrand-canonical ensemble. Since acceptance probabilities for direct transitions between DMPC and DHPC were negligible, a third component with intermediate tail length (didecanoylphosphatidylcholine, di-C10 tails) was included at a low concentration to serve as an intermediate for transitions between DMPC and DHPC. Strong enrichment of DHPC is seen at ribbon and pore edges, with an excess linear density of ∼3 nm−1. The simulation model yields estimates for the onset of edge stability with increasing bilayer DHPC content between 5% and 15% DHPC at 300 K and between 7% and 17% DHPC at 323 K, higher than experimental estimates. Local structure and composition at points of close contact between pores suggest a possible mechanism for effective attractions between pores, providing a rationalization for the tendency of bicelle mixtures to aggregate into perforated vesicles and perforated sheets.

The partitioning of lipids among different microenvironments in a bilayer is of considerable relevance to characterization of composition variations in biomembranes. Atomistic simulation has been ill-suited to model equilibrated lipid mixtures because the time required for diffusive exchange of lipids among microenvironments exceeds typical submicrosecond molecular dynamics trajectories. A method to facilitate local composition fluctuations, using Monte Carlo mutations to change lipid structures within the semigrand-canonical ensemble (at a fixed difference in component chemical potentials, Δμ), was recently implemented to address this challenge. This technique was applied here to mixtures of dimyristoylphosphatidylcholine and a shorter-tail lipid, either symmetric (didecanoylphosphatidylcholine (DDPC)) or asymmetric (hexanoyl-myristoylphosphatidylcholine), arranged in two types of structure: bilayer ribbons and buckled bilayers. In ribbons, the shorter-tail component showed a clear enrichment at the highly curved rim, more so for hexanoyl-myristoylphosphatidylcholine than for DDPC. Results on buckled bilayers were variable. Overall, the DDPC content of buckled bilayers tended to exceed by several percent the DDPC content of flat ones simulated at the same Δμ, but only for mixtures with low overall DDPC content. Within the buckled bilayer structure, no correlation could be resolved between the sign or magnitude of the local curvature of a leaflet and the mean local lipid composition. Results are discussed in terms of packing constraints, surface area/volume ratios, and curvature elasticity.