Bringing discrete-time Langevin splitting methods into agreement with thermodynamics.

In light of the recently published complete set of statistically correct Grønbech-Jensen (GJ) methods for discrete-time thermodynamics, we revise a differential operator splitting method for the Langevin equation in order to comply with the basic GJ thermodynamic sampling features, namely, the Boltzmann distribution and Einstein diffusion, in linear systems. This revision, which is based on the introduction of time scaling along with flexibility of a discrete-time velocity attenuation parameter, provides a direct link between the ABO splitting formalism and the GJ methods. This link brings about the conclusion that any GJ method has at least weak second order accuracy in the applied time step. It further helps identify a novel half-step velocity, which simultaneously produces both correct kinetic statistics and correct transport measures for any of the statistically sound GJ methods. Explicit algorithmic expressions are given for the integration of the new half-step velocity into the GJ set of methods. Numerical simulations, including quantum-based molecular dynamics (QMD) using the QMD suite Los Alamos Transferable Tight-Binding for Energetics, highlight the discussed properties of the algorithms as well as exhibit the direct application of robust, time-step-independent stochastic integrators to QMD.

Coexistence of Lipid Phases Stabilizes Interstitial Water in the Outer Layer of Mammalian Skin.

The lipid matrix in the outer layer of mammalian skin, the stratum corneum, has been previously investigated by multiple biophysical techniques aimed at identifying hydrophilic and lipophilic pathways of permeation. Although consensus is developing over the microscopic structure of the lipid matrix, no molecular-resolution model describes the permeability of all chemical species simultaneously. Using molecular dynamics simulations of a model mixture of skin lipids, the self-assembly of the lipid matrix lamellae has been studied. At higher humidity, the resulting lamellar phase is maintained by partitioning excess water into isolated droplets of controlled size and spatial distribution. The droplets may fuse together to form intralamellar water channels, thereby providing a pathway for the permeation of hydrophilic species. These results reconcile competing data on the outer skin's structure and broaden the scope of molecular-based methods to improve the safety of topical products and to advance transdermal drug delivery.

Direct Derivation of Free Energies of Membrane Deformation and Other Solvent Density Variations From Enhanced Sampling Molecular Dynamics.

We report a methodology to calculate the free energy of a shape transformation in a lipid membrane directly from a molecular dynamics simulation. The bilayer need not be homogeneous or symmetric and can be atomically detailed or coarse grained. The method is based on a collective variable that quantifies the similarity between the membrane and a set of predefined density distributions. Enhanced sampling of this "Multi-Map" variable re-shapes the bilayer and permits the derivation of the corresponding potential of mean force. Calculated energies thus reflect the dynamic interplay of atoms and molecules, rather than postulated effects. Evaluation of deformations of different shape, amplitude, and range demonstrates that the macroscopic bending modulus assumed by the Helfrich-Canham model is increasingly unsuitable below the 100-Å scale. In this range of major biological significance, direct free-energy calculations reveal a much greater plasticity. We also quantify the stiffening effect of cholesterol on bilayers of different composition and compare with experiments. Lastly, we illustrate how this approach facilitates analysis of other solvent reorganization processes, such as hydrophobic hydration. Published 2019. This article is a U.S. Government work and is in the public domain in the USA.

The challenge of stochastic Størmer-Verlet thermostats generating correct statistics.

In light of the recently developed complete GJ set of single random variable stochastic, discrete-time Størmer-Verlet algorithms for statistically accurate simulations of Langevin equations [N. Grønbech-Jensen, Mol. Phys. 118, e1662506 (2020)], we investigate two outstanding questions: (1) Are there any algorithmic or statistical benefits from including multiple random variables per time step and (2) are there objective reasons for using one or more methods from the available set of statistically correct algorithms? To address the first question, we assume a general form for the discrete-time equations with two random variables and then follow the systematic, brute-force GJ methodology by enforcing correct thermodynamics in linear systems. It is concluded that correct configurational Boltzmann sampling of a particle in a harmonic potential implies correct configurational free-particle diffusion and that these requirements only can be accomplished if the two random variables per time step are identical. We consequently submit that the GJ set represents all possible stochastic Størmer-Verlet methods that can reproduce time step-independent statistics of linear systems. The second question is thus addressed within the GJ set. Based on numerical simulations of complex molecular systems, as well as on analytic considerations, we analyze apparent friction-induced differences in the stability of the methods. We attribute these differences to an inherent, friction-dependent discrete-time scaling, which depends on the specific method. We suggest that the method with the simplest interpretation of temporal scaling, the GJ-I/GJF-2GJ method, be preferred for statistical applications.

Data-guided Multi-Map variables for ensemble refinement of molecular movies.

Driving molecular dynamics simulations with data-guided collective variables offer a promising strategy to recover thermodynamic information from structure-centric experiments. Here, the three-dimensional electron density of a protein, as it would be determined by cryo-EM or x-ray crystallography, is used to achieve simultaneously free-energy costs of conformational transitions and refined atomic structures. Unlike previous density-driven molecular dynamics methodologies that determine only the best map-model fits, our work employs the recently developed Multi-Map methodology to monitor concerted movements within equilibrium, non-equilibrium, and enhanced sampling simulations. Construction of all-atom ensembles along the chosen values of the Multi-Map variable enables simultaneous estimation of average properties, as well as real-space refinement of the structures contributing to such averages. Using three proteins of increasing size, we demonstrate that biased simulation along the reaction coordinates derived from electron densities can capture conformational transitions between known intermediates. The simulated pathways appear reversible with minimal hysteresis and require only low-resolution density information to guide the transition. The induced transitions also produce estimates for free energy differences that can be directly compared to experimental observables and population distributions. The refined model quality is superior compared to those found in the Protein Data Bank. We find that the best quantitative agreement with experimental free-energy differences is obtained using medium resolution density information coupled to comparatively large structural transitions. Practical considerations for probing the transitions between multiple intermediate density states are also discussed.

Scalable molecular dynamics on CPU and GPU architectures with NAMD.

NAMDis a molecular dynamics program designed for high-performance simulations of very large biological objects on CPU- and GPU-based architectures. NAMD offers scalable performance on petascale parallel supercomputers consisting of hundreds of thousands of cores, as well as on inexpensive commodity clusters commonly found in academic environments. It is written in C++ and leans on Charm++ parallel objects for optimal performance on low-latency architectures. NAMD is a versatile, multipurpose code that gathers state-of-the-art algorithms to carry out simulations in apt thermodynamic ensembles, using the widely popular CHARMM, AMBER, OPLS, and GROMOS biomolecular force fields. Here, we review the main features of NAMD that allow both equilibrium and enhanced-sampling molecular dynamics simulations with numerical efficiency. We describe the underlying concepts utilized by NAMD and their implementation, most notably for handling long-range electrostatics; controlling the temperature, pressure, and pH; applying external potentials on tailored grids; leveraging massively parallel resources in multiple-copy simulations; and hybrid quantum-mechanical/molecular-mechanical descriptions. We detail the variety of options offered by NAMD for enhanced-sampling simulations aimed at determining free-energy differences of either alchemical or geometrical transformations and outline their applicability to specific problems. Last, we discuss the roadmap for the development of NAMD and our current efforts toward achieving optimal performance on GPU-based architectures, for pushing back the limitations that have prevented biologically realistic billion-atom objects to be fruitfully simulated, and for making large-scale simulations less expensive and easier to set up, run, and analyze. NAMD is distributed free of charge with its source code at www.ks.uiuc.edu.

High-resolution structures of the M2 channel from influenza A virus reveal dynamic pathways for proton stabilization and transduction.

The matrix 2 (M2) protein from influenza A virus is a proton channel that uses His37 as a selectivity filter. Here we report high-resolution (1.10 Å) cryogenic crystallographic structures of the transmembrane domain of M2 at low and high pH. These structures reveal that waters within the pore form hydrogen-bonded networks or "water wires" spanning 17 Å from the channel entrance to His37. Pore-lining carbonyl groups are well situated to stabilize hydronium via second-shell interactions involving bridging water molecules. In addition, room temperature crystallographic structures indicate that water becomes increasingly fluid with increasing temperature and decreasing pH, despite the higher electrostatic field. Complementary molecular dynamics simulations reveal a collective switch of hydrogen bond orientations that can contribute to the directionality of proton flux as His37 is dynamically protonated and deprotonated in the conduction cycle.

Structural Refinement of Proteins by Restrained Molecular Dynamics Simulations with Non-interacting Molecular Fragments.

The knowledge of multiple conformational states is a prerequisite to understand the function of membrane transport proteins. Unfortunately, the determination of detailed atomic structures for all these functionally important conformational states with conventional high-resolution approaches is often difficult and unsuccessful. In some cases, biophysical and biochemical approaches can provide important complementary structural information that can be exploited with the help of advanced computational methods to derive structural models of specific conformational states. In particular, functional and spectroscopic measurements in combination with site-directed mutations constitute one important source of information to obtain these mixed-resolution structural models. A very common problem with this strategy, however, is the difficulty to simultaneously integrate all the information from multiple independent experiments involving different mutations or chemical labels to derive a unique structural model consistent with the data. To resolve this issue, a novel restrained molecular dynamics structural refinement method is developed to simultaneously incorporate multiple experimentally determined constraints (e.g., engineered metal bridges or spin-labels), each treated as an individual molecular fragment with all atomic details. The internal structure of each of the molecular fragments is treated realistically, while there is no interaction between different molecular fragments to avoid unphysical steric clashes. The information from all the molecular fragments is exploited simultaneously to constrain the backbone to refine a three-dimensional model of the conformational state of the protein. The method is illustrated by refining the structure of the voltage-sensing domain (VSD) of the Kv1.2 potassium channel in the resting state and by exploring the distance histograms between spin-labels attached to T4 lysozyme. The resulting VSD structures are in good agreement with the consensus model of the resting state VSD and the spin-spin distance histograms from ESR/DEER experiments on T4 lysozyme are accurately reproduced.

Infrared and fluorescence assessment of the hydration status of the tryptophan gate in the influenza A M2 proton channel.

The M2 proton channel of the influenza A virus has been the subject of extensive studies because of its critical role in viral replication. As such, we now know a great deal about its mechanism of action, especially how it selects and conducts protons in an asymmetric fashion. The conductance of this channel is tuned to conduct protons at a relatively low biologically useful rate, which allows acidification of the viral interior of a virus entrapped within an endosome, but not so great as to cause toxicity to the infected host cell prior to packaging of the virus. The dynamic, structural and chemical features that give rise to this tuning are not fully understood. Herein, we use a tryptophan (Trp) analog, 5-cyanotryptophan, and various methods, including linear and nonlinear infrared spectroscopies, static and time-resolved fluorescence techniques, and molecular dynamics simulations, to site-specifically interrogate the structure and hydration dynamics of the Trp41 gate in the transmembrane domain of the M2 proton channel. Our results suggest that the Trp41 sidechain adopts the t90 rotamer, the χ2 dihedral angle of which undergoes an increase of approximately 35° upon changing the pH from 7.4 to 5.0. Furthermore, we find that Trp41 is situated in an environment lacking bulk-like water, and somewhat surprisingly, the water density and dynamics do not show a measurable difference between the high (7.4) and low (5.0) pH states. Since previous studies have shown that upon channel opening water flows into the cavity above the histidine tetrad (His37), the present finding thus provides evidence indicating that the lack of sufficient water molecules near Trp41 needed to establish a continuous hydrogen bonding network poses an additional energetic bottleneck for proton conduction.

Pore waters regulate ion permeation in a calcium release-activated calcium channel.

The recent crystal structure of Orai, the pore unit of a calcium release-activated calcium (CRAC) channel, is used as the starting point for molecular dynamics and free-energy calculations designed to probe this channel's conduction properties. In free molecular dynamics simulations, cations localize preferentially at the extracellular channel entrance near the ring of Glu residues identified in the crystal structure, whereas anions localize in the basic intracellular half of the pore. To begin to understand ion permeation, the potential of mean force (PMF) was calculated for displacing a single Na(+) ion along the pore of the CRAC channel. The computed PMF indicates that the central hydrophobic region provides the major hindrance for ion diffusion along the permeation pathway, thereby illustrating the nonconducting nature of the crystal structure conformation. Strikingly, further PMF calculations demonstrate that the mutation V174A decreases the free energy barrier for conduction, rendering the channel effectively open. This seemingly dramatic effect of mutating a nonpolar residue for a smaller nonpolar residue in the pore hydrophobic region suggests an important role for the latter in conduction. Indeed, our computations show that even without significant channel-gating motions, a subtle change in the number of pore waters is sufficient to reshape the local electrostatic field and modulate the energetics of conduction, a result that rationalizes recent experimental findings. The present work suggests the activation mechanism for the wild-type CRAC channel is likely regulated by the number of pore waters and hence pore hydration governs the conductance.

Structure and inhibition of the drug-resistant S31N mutant of the M2 ion channel of influenza A virus.

The influenza A virus M2 proton channel (A/M2) is the target of the antiviral drugs amantadine and rimantadine, whose use has been discontinued due to widespread drug resistance. Among the handful of drug-resistant mutants, S31N is found in more than 95% of the currently circulating viruses and shows greatly decreased inhibition by amantadine. The discovery of inhibitors of S31N has been hampered by the limited size, polarity, and dynamic nature of its amantadine-binding site. Nevertheless, we have discovered small-molecule drugs that inhibit S31N with potencies greater than amantadine's potency against WT M2. Drug binding locks the protein into a well-defined conformation, and the NMR structure of the complex shows the drug bound in the homotetrameric channel, threaded between the side chains of Asn31. Unrestrained molecular dynamics simulations predicted the same binding site. This S31N inhibitor, like other potent M2 inhibitors, contains a charged ammonium group. The ammonium binds as a hydrate to one of three sites aligned along the central cavity that appear to be hotspots for inhibition. These sites might stabilize hydronium-like species formed as protons diffuse through the outer channel to the proton-shuttling residue His37 near the cytoplasmic end of the channel.

How Sensitive is the Amide†I Vibration of the Polypeptide Backbone to Electric Fields?

Site-selective isotopic labeling of amide carbonyls offers a nonperturbative means to introduce a localized infrared probe into proteins. Although this strategy has been widely used to investigate various biological questions, the dependence of the underlying amide I vibrational frequency on electric fields (or Stark tuning rate) has not been fully determined, which prevents it from being used in a quantitative manner in certain applications. Herein, through the use of experiments and molecular dynamics simulations, the Stark tuning rate of the amide I vibration of an isotopically labeled backbone carbonyl in a transmembrane α-helix is determined to be approximately 1.4 cm(-1) /(MV/cm). This result provides a quantitative basis for using this vibrational model to assess local electric fields in proteins, among other applications. For instance, by using this value, we are able to show that the backbone region of a dipeptide has a surprisingly low dielectric constant.

Molecular dynamics simulations of cholesterol-rich membranes using a coarse-grained force field for cyclic alkanes.

The architecture of a biological membrane hinges upon the fundamental fact that its properties are determined by more than the sum of its individual components. Studies on model membranes have shown the need to characterize in molecular detail how properties such as thickness, fluidity, and macroscopic bending rigidity are regulated by the interactions between individual molecules in a non-trivial fashion. Simulation-based approaches are invaluable to this purpose but are typically limited to short sampling times and model systems that are often smaller than the required properties. To alleviate both limitations, the use of coarse-grained (CG) models is nowadays an established computational strategy. We here present a new CG force field for cholesterol, which was developed by using measured properties of small molecules, and can be used in combination with our previously developed force field for phospholipids. The new model performs with precision comparable to atomistic force fields in predicting the properties of cholesterol-rich phospholipid bilayers, including area per lipid, bilayer thickness, tail order parameter, increase in bending rigidity, and propensity to form liquid-ordered domains in ternary mixtures. We suggest the use of this model to quantify the impact of cholesterol on macroscopic properties and on microscopic phenomena involving localization and trafficking of lipids and proteins on cellular membranes.

Flipping in the pore: discovery of dual inhibitors that bind in different orientations to the wild-type versus the amantadine-resistant S31N mutant of the influenza A virus M2 proton channel.

Influenza virus infections lead to numerous deaths and millions of hospitalizations each year. One challenge facing anti-influenza drug development is the heterogeneity of the circulating influenza viruses, which comprise several strains with variable susceptibility to antiviral drugs. For example, the wild-type (WT) influenza A viruses, such as the seasonal H1N1, tend to be sensitive to antiviral drugs, amantadine and rimantadine, while the S31N mutant viruses, such as the pandemic 2009 H1N1 (H1N1pdm09) and seasonal H3N2, are resistant to this class of drugs. Thus, drugs targeting both WT and the S31N mutant are highly desired. We report our design of a novel class of dual inhibitors along with their ion channel blockage and antiviral activities. The potency of the most active compound 11 in inhibiting WT and the S31N mutant influenza viruses is comparable with that of amantadine in inhibiting WT influenza virus. Solution NMR studies and molecular dynamics (MD) simulations of drug-M2 interactions supported our design hypothesis: namely, the dual inhibitor binds in the WT M2 channel with an aromatic group facing down toward the C-terminus, while the same drug binds in the S31N M2 channel with its aromatic group facing up toward the N-terminus. The flip-flop mode of drug binding correlates with the structure-activity relationship (SAR) and has paved the way for the next round of rational design of broad-spectrum antiviral drugs.

Structure of water at charged interfaces: a molecular dynamics study.

The properties of water molecules located close to an interface deviate significantly from those observed in the homogeneous bulk liquid. The length scale over which this structural perturbation persists (the so-called interfacial depth) is the object of extensive investigations. The situation is particularly complicated in the presence of surface charges that can induce long-range orientational ordering of water molecules, which in turn dictate diverse processes, such as mineral dissolution, heterogeneous catalysis, and membrane chemistry. To characterize the fundamental properties of interfacial water, we performed molecular dynamics (MD) simulations on alkali chloride solutions in the presence of two types of idealized charged surfaces: one with the charge density localized at discrete sites and the other with a homogeneously distributed charge density. We find that, in addition to a diffuse region where water orientation shows no layering, the interface region consists of a "compact layer" of solvent next to the surface that is not described in classical electric double layer theories. The depth of the diffuse solvent layer is sensitive to the type of charge distributions on the surface and the ionic strength. Simulations of the aqueous interface of a realistic model of negatively charged amorphous silica show that the water orientation and the distribution of ions strongly depend on the identity of the cations (Na(+) vs Cs(+)) and are not well represented by a simplistic homogeneous charge distribution model. While the compact layer shows different solvent net orientation and depth for Na(+) vs Cs(+), the depth (~1 nm) of the diffuse layer of oriented waters is independent of the identity of the cation screening the charge. The details of interfacial water orientation revealed here go beyond the traditionally used double and triple layer models and provide a microscopic picture of the aqueous/mineral interface that complements recent surface specific experimental studies.

Dehydration of multilamellar fatty acid membranes: towards a computational model of the stratum corneum.

The level of hydration controls the cohesion between apposed lamellae of saturated free fatty acids found in the lipid matrix of stratum corneum, the outermost layer of mammalian skin. This multilamellar lipid matrix is highly impermeable to water and ions, so that the local hydration shell of its fatty acids may not always be in equilibrium with the acidity and relative humidity, which significantly change over a course of days during skin growth. The homeostasis of the stratum corneum at each moment of its growth likely requires a balance between two factors, which affect in opposite ways the diffusion of hydrophilic species through the stratum corneum: (i) an increase in water order as the lipid lamellae come in closer contact, and (ii) a decrease in water order as the fraction of charged fatty acids is lowered by pH. Herein molecular dynamics simulations are employed to estimate the impact of both effects on water molecules confined between lamellae of fatty acids. Under conditions where membrane undulations are energetically favorable, the charged fatty acids are able to sequester cations around points of contact between lamellae that are fully dehydrated, while essentially maintaining a multilamellar structure for the entire system. This observation suggests that the undulations of the fatty acid lamellae control the diffusion of hydrophilic species through the water phase by altering the positional and rotational order of water molecules in the embedded/occluded "droplets."

Transmembrane orientation and possible role of the fusogenic peptide from parainfluenza virus 5 (PIV5) in promoting fusion.

Membrane fusion is required for diverse biological functions ranging from viral infection to neurotransmitter release. Fusogenic proteins increase the intrinsically slow rate of fusion by coupling energetically downhill conformational changes of the protein to kinetically unfavorable fusion of the membrane-phospholipid bilayers. Class I viral fusogenic proteins have an N-terminal hydrophobic fusion peptide (FP) domain, important for interaction with the target membrane, plus a C-terminal transmembrane (C-term-TM) helical membrane anchor. The role of the water-soluble regions of fusogenic proteins has been extensively studied, but the contributions of the membrane-interacting FP and C-term-TM peptides are less well characterized. Typically, FPs are thought to bind to membranes at an angle that allows helix penetration but not traversal of the lipid bilayer. Here, we show that the FP from the paramyxovirus parainfluenza virus 5 fusogenic protein, F, forms an N-terminal TM helix, which self-associates into a hexameric bundle. This FP also interacts strongly with the C-term-TM helix. Thus, the fusogenic F protein resembles SNARE proteins involved in vesicle fusion by having water-soluble coiled coils that zipper during fusion and TM helices in both membranes. By analogy to mechanosensitive channels, the force associated with zippering of the water-soluble coiled-coil domain is expected to lead to tilting of the FP helices, promoting interaction with the C-term-TM helices. The energetically unfavorable dehydration of lipid headgroups of opposing bilayers is compensated by thermodynamically favorable interactions between the FP and C-term-TM helices as the coiled coils zipper into the membrane phase, leading to a pore lined by both lipid and protein.

Structure and mechanism of proton transport through the transmembrane tetrameric M2 protein bundle of the influenza A virus.

The M2 proton channel from influenza A virus is an essential protein that mediates transport of protons across the viral envelope. This protein has a single transmembrane helix, which tetramerizes into the active channel. At the heart of the conduction mechanism is the exchange of protons between the His37 imidazole moieties of M2 and waters confined to the M2 bundle interior. Protons are conducted as the total charge of the four His37 side chains passes through 2(+) and 3(+) with a pK(a) near 6. A 1.65 A resolution X-ray structure of the transmembrane protein (residues 25-46), crystallized at pH 6.5, reveals a pore that is lined by alternating layers of sidechains and well-ordered water clusters, which offer a pathway for proton conduction. The His37 residues form a box-like structure, bounded on either side by water clusters with well-ordered oxygen atoms at close distance. The conformation of the protein, which is intermediate between structures previously solved at higher and lower pH, suggests a mechanism by which conformational changes might facilitate asymmetric diffusion through the channel in the presence of a proton gradient. Moreover, protons diffusing through the channel need not be localized to a single His37 imidazole, but instead may be delocalized over the entire His-box and associated water clusters. Thus, the new crystal structure provides a possible unification of the discrete site versus continuum conduction models.

Membrane free-energy landscapes derived from atomistic dynamics explain nonuniversal cholesterol-induced stiffening.

All lipid membranes have inherent morphological preferences and resist deformation. Yet adaptations in membrane shape can and do occur at multiple length scales. While this plasticity is crucial for cellular physiology, the factors controlling the morphological energetics of lipid bilayers and the dominant mechanisms of membrane remodeling remain unclear. An ongoing debate regarding the universality of the stiffening effect of cholesterol underscores the challenges facing this field, both experimentally and theoretically, even for simple lipid mixtures. On the computational side, we have argued that enhanced- sampling all-atom molecular dynamics simulations are uniquely suited for quantification of membrane conformational energetics, not only because they minimize a-priori assumptions, but also because they permit analysis of bilayers in deformed states. To showcase this approach, we examine reported inconsistencies between alternative experimental measurements of bending moduli for cholesterol-enriched membranes. Specifically, we analyze lipid bilayers with different chain saturation, and compute free-energy landscapes for curvature deformations distributed over areas from ∼5 to ∼60 nm 2 . These enhanced simulations, totaling over 100 microseconds of sampling time, enable us to directly quantify both bending and tilt moduli, and to dissect the contributing factors and molecular mechanisms of curvature generation at each length scale. Our results show that cholesterol effects are lipid-specific, in agreement with giantvesicle measurements, and explain why experiments probing nanometer scale lipid dynamics diverge. In summary, we demonstrate that quantitative structure-mechanics relationships can now be established for heterogenous membranes, paving the way for addressing open fundamental questions in cell membrane mechanics. Significance: Elucidating the energetics and mechanisms of membrane remodeling is an essential step towards understanding cell physiology. This problem is challenging, however, because membrane bending involves both large-scale and atomic-level dynamics, which are difficult to measure simultaneously. A recent controversy regarding the stiffening effect of cholesterol, which is ubiquitous in animal cells, illustrates this challenge. We show how enhanced molecular-dynamics simulations can bridge this length-scale gap and reconcile seemingly incongruent observations. This approach facilitates a conceptual connection between lipid chemistry and membrane mechanics, thereby providing a solid basis for future research on remodeling phenomena, such as in membrane trafficking or viral infection.

Membrane free-energy landscapes derived from atomistic dynamics explain nonuniversal cholesterol-induced stiffening.

All lipid membranes have inherent morphological preferences and resist deformation. Yet adaptations in membrane shape can and do occur at multiple length scales. While this plasticity is crucial for cellular physiology, the factors controlling the morphological energetics of lipid bilayers and the dominant mechanisms of membrane remodeling remain to be fully understood. An ongoing debate regarding the universality of the stiffening effect of cholesterol underscores the challenges facing this field, both experimentally and theoretically, even for simple lipid mixtures. On the computational side, we have argued that enhanced-sampling all-atom molecular dynamics simulations are uniquely suited for the quantification of membrane conformational energetics, as they minimize a priori assumptions and permit analysis of bilayers in deformed states. To showcase this approach, we examine reported inconsistencies between alternative experimental measurements of bending moduli for cholesterol-enriched membranes. Specifically, we analyze lipid bilayers with different chain saturation and compute free-energy landscapes for curvature deformations distributed over areas from ∼5 to ∼60 nm2. These enhanced simulations, totaling over 100 µs of sampling time, enable us to directly quantify both bending and tilt moduli and to dissect the contributing factors and molecular mechanisms of curvature generation at each length scale. Our results show that the effects of cholesterol on bending rigidity are lipid-specific and suggest that this specificity arises from differences in the torsional dynamics of the acyl chains. In summary, we demonstrate that quantitative relationships can now be established between lipid structure and bending energetics, paving the way for addressing open fundamental questions in cell membrane mechanics.