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.

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.

SPICA Force Field for Proteins and Peptides

A coarse-grained (CG) model for peptides and proteins was developed as an extension of the Surface Property fItting Coarse grAined (SPICA) force field (FF). The model was designed to examine membrane proteins that are fully compatible with the lipid membranes of the SPICA FF. A preliminary version of this protein model was created using thermodynamic properties, including the surface tension and density in the SPICA (formerly called SDK) FF. In this study, we improved the CG protein model to facilitate molecular dynamics (MD) simulations with a reproduction of multiple properties from both experiments and all-atom (AA) simulations. An elastic network model was adopted to maintain the secondary structure within a single chain. The side-chain analogues reproduced the transfer free energy profiles across the lipid membrane and demonstrated reasonable association free energy (potential of mean force) in water compared to those from AA MD. A series of peptides/proteins adsorbed onto or penetrated into the membrane simulated by the CG MD correctly predicted the penetration depths and tilt angles of peripheral and transmembrane peptides/proteins as comparable to those in the orientations of proteins in membranes (OPM) database. In addition, the dimerization free energies of several transmembrane helices within a lipid bilayer were comparable to those from experimental estimation. Application studies on a series of membrane protein assemblies, scramblases, and poliovirus capsids demonstrated the good performance of the SPICA FF.

Human Learning for Molecular Simulations: The Collective Variables Dashboard in VMD

The Collective Variables Dashboard is a software tool for real-time, seamless exploration of molecular structures and trajectories in a customizable space of collective variables. The Dashboard arises from the integration of the Collective Variables Module (also known as Colvars) with the visualization software VMD, augmented with a fully discoverable graphical interface offering interactive workflows for the design and analysis of collective variables. Typical use cases include a priori design of collective variables for enhanced sampling and free energy simulations as well as analysis of any type of simulation or collection of structures in a collective variable space. A combination of those cases commonly occurs when preliminary simulations, biased or unbiased, reveal that an optimized set of collective variables is necessary to improve sampling in further simulations. Then the Dashboard provides an efficient way to intuitively explore the space of likely collective variables, validate them on existing data, and use the resulting collective variable definitions directly in further biased simulations using the Collective Variables Module. Visualization of biasing energies and forces is proposed to help analyze or plan biased simulations. We illustrate the use of the Dashboard on two applications: discovering coordinates to describe ligand unbinding from a protein binding site and designing volume-based variables to bias the hydration of a transmembrane pore.

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.

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.

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.

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.

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.

Large-scale state-dependent membrane remodeling by a transporter protein

That channels and transporters can influence the membrane morphology is increasingly recognized. Less appreciated is that the extent and free-energy cost of these deformations likely varies among different functional states of a protein, and thus, that they might contribute significantly to defining its mechanism. We consider the trimeric Na⁺-aspartate symporter Glt_Ph, a homolog of an important class of neurotransmitter transporters, whose mechanism entails one of the most drastic structural changes known. Molecular simulations indicate that when the protomers become inward-facing, they cause deep, long-ranged, and yet mutually-independent membrane deformations. Using a novel simulation methodology, we estimate that the free-energy cost of this membrane perturbation is in the order of 6–7 kcal/mol per protomer. Compensating free-energy contributions within the protein or its environment must thus stabilize this inward-facing conformation for the transporter to function. We discuss these striking results in the context of existing experimental observations for this and other transporters.

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.

Structural Characterization of Biomolecules through Atomistic Simulations Guided by DEER Measurements

Double electron-electron resonance (DEER) is a popular technique that exploits attached spin labels to probe the collective dynamics of biomolecules in a native environment. Like most spectroscopic approaches, DEER detects an ensemble of states accounting for biomolecular dynamics as well as the labels' intrinsic flexibility. Hence, the DEER data alone do not provide high-resolution structural information. To disentangle this variability, we introduce a minimally biased simulation method to sample a structural ensemble that reproduces multiple experimental signals within the uncertainty. In contrast to previous approaches, our method targets the raw data themselves, and thereby it brings forth an unbiased molecular interpretation of the experiments. After validation on the T4 lysozyme, we apply this technique to interpret recent DEER experiments on a membrane transporter binding protein (VcSiaP). The results highlight the large-scale conformational movement that occurs upon substrate binding and reveal that the unbound VcSiaP is more open in solution than the X-ray structure.

Mechanically Strong Polymer Sheets from Aligned Ultrahigh-Molecular-Weight Polyethylene Nanocomposites

Ultrahigh-molecular-weight polyethylene (UHMWPE) is of great interest as a next-generation body armor material because of its superior mechanical properties. However, such unique properties depend critically on its microscopic structure characteristics, including the degree of crystallinity, chain alignment, and morphology. Here, we present a highly aligned UHMWPE and its composite sheets containing uniformly dispersed boron nitride (BN) nanosheets. The dispersion of BN nanosheets into the UHMWPE matrix increases its mechanical properties over a broad temperature range. Experiments and simulation confirm that the alignment of chain segments in the composite matrix increases with temperature, leading to an improvement in mechanical properties at high temperature. Together with the large thermal conductivity of UHMWPE and BN, our findings serve to expand the application spectrum of highly aligned polymer nanocomposite materials for ballistic panels and body armor over a broad range of temperatures.

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.

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.

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.

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."

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.

Hydrogen-bonded water molecules in the M2 channel of the influenza A virus guide the binding preferences of ammonium-based inhibitors

The tetrameric M2 proton channel of influenza A virus is an integral membrane protein responsible for the acidification of the viral interior. Drugs such as amantadine target the transmembrane region of wild type M2 by acting as pore blockers. However, a number of mutations affecting this domain confer drug resistance, prompting the need for alternative inhibitors. The availability of high-resolution structures of drug-bound M2, paired with computational investigations, revealed that inhibitors can bind at different sites, and provided useful insights in understanding the principles governing proton conduction. Here, we investigated by computation the energetic and geometric factors determining the relative stability of pore blockers at individual sites of different M2 strains. We found that local free energy minima along the translocation pathway of positively charged chemical species correspond to experimentally determined binding sites of inhibitors. Then, by examining the structure of water clusters hydrating each site, as well as of those displaced by binding of hydrophobic scaffolds, we predicted the binding preferences of M2 ligands. This information can be used to guide the identification of novel classes of inhibitors.

Proton release from the histidine-tetrad in the M2 channel of the influenza A virus

The activity of the M2 proton channel of the influenza A virus is controlled by pH. The tautomeric state and conformation of His37, a key residue in the M2 transmembrane four-helix bundle, controls the gating of the channel. Previously, we compared the energetics and dynamics of two alternative conformations of the doubly protonated state at neutral pH, namely, a 4-fold symmetric “histidine-box” and a 2-fold symmetric “dimer-of-dimers”, and proposed a multiconfiguration model for this charge state. Here, we elaborate this model by further studying configurations of the His37 tetrad in the triply protonated state and its subsequent deprotonation via quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations, starting with the aforementioned configurations, to gain information about proton release in a viral membrane environment. Interestingly, the two configurations converge under acidic pH conditions. Protons can be transferred from one charged His37 to a neighboring water cluster at the C-terminal side of the channel when the Trp41 gate is open transiently. With limited backbone expansion, the free energy barrier for proton release to the viral interior at low pH is ∼6.5 kcal/mol in both models, which is much lower than at either neutral pH or for an isolated His37 cluster without a membrane environment. Our calculations also suggest that the M2 protein would seem to exclude the entrance of anions into the central channel through a special mechanism, due to the latter’s potential inhibitory effect on proton conduction.

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.

Exploring Histidine Conformations in the M2 Channel Lumen of the Influenza A Virus at Neutral pH via Molecular Simulations

The pH-regulated M2 proton channel from the influenza A virus has a His-tetrad in its transmembrane (TM) domain that is essential for proton conduction. At neutral pH, the tetrad has been observed in two distinct configurations, the "His-box" and "dimer-of-dimers", with similar backbone structures suggesting competing models for proton conduction. Here, we propose that both conformations can play a role. In support of this hypothesis, we used molecular dynamics simulations based on density functional theory to simulate the M2-TM domain and force-field-based simulations to estimate the relevant free-energy barriers. Both configurations are stable on accessible simulation time scales, and transitions between them occur faster than the millisecond time scale of proton conduction. Moreover, the deprotonation energy is too high for spontaneous conduction, consistent with their occurrence in the low-current regime. Our computations support a multiconfiguration model with different population levels, thereby connecting experimental data obtained under different conditions.

Discovery of novel dual inhibitors of the wild-type and the most prevalent drug-resistant mutant, S31N, of the M2 proton channel from influenza A virus

Anti-influenza drugs, amantadine and rimantadine, targeting the M2 channel from influenza A virus are no longer effective because of widespread drug resistance. S31N is the predominant and amantadine-resistant M2 mutant, present in almost all of the circulating influenza A strains as well as in the pandemic 2009 H1N1 and the highly pathogenic H5N1 flu strains. Thus, there is an urgent need to develop second-generation M2 inhibitors targeting the S31N mutant. However, the S31N mutant presents a huge challenge to drug discovery, and it has been considered undruggable for several decades. Using structural information, classical medicinal chemistry approaches, and M2-specific biological testing, we discovered benzyl-substituted amantadine derivatives with activity against both S31N and WT, among which 4-(adamantan-1-ylaminomethyl)-benzene-1,3-diol (44) is the most potent dual inhibitor. These inhibitors demonstrate that S31N is a druggable target and provide a new starting point to design novel M2 inhibitors that address the problem of drug-resistant influenza A infections.

Molecular dynamics simulation directed rational design of inhibitors targeting drug-resistant mutants of influenza A virus M2

Influenza A virus M2 (A/M2) forms a homotetrameric proton selective channel in the viral membrane. It has been the drug target of antiviral drugs such as amantadine and rimantadine. However, most of the current virulent influenza A viruses carry drug-resistant mutations alongside the drug binding site, such as S31N, V27A, and L26F, etc., each of which might be dominant in a given flu season. Among these mutations, the V27A mutation was prevalent among transmissible viruses under drug selection pressure. Until now, V27A has not been successfully targeted by small molecule inhibitors, despite years of extensive medicinal chemistry research efforts and high throughput screening. Guided by molecular dynamics (MD) simulation of drug binding and the influence of drug binding on the dynamics of A/M2 from earlier experimental studies, we designed a series of potent spirane amine inhibitors targeting not only WT, but also both A/M2-27A and L26F mutants with IC(50)s similar to that seen for amantadine's inhibition of the WT channel. The potencies of these inhibitors were further demonstrated in experimental binding and plaque reduction assays. These results demonstrate the power of MD simulations to probe the mechanism of drug binding as well as the ability to guide design of inhibitors of targets that had previously appeared to be undruggable.

Structural biology. The flu's proton escort

A flurry of structural data provides sometimes conflicting insights into the M2 proton channel.

Multiple Proton Confinement in the M2 Channel from the Influenza A Virus

The tetrameric M2 protein bundle of the influenza A virus is the proton channel responsible for the acidification of the viral interior, a key step in the infection cycle. Selective proton transport is achieved by successive protonation of the conserved histidine amino acids at position 37. A recent X-ray structure of the tetrameric transmembrane (TM) domain of the protein (residues 22-46) resolved several water clusters in the channel lumen, which suggest possible proton pathways to the His37 residues. To explore this hypothesis, we have carried out molecular dynamics (MD) simulations of a proton traveling towards the His37 side chains using MD with classical and quantum force fields. Diffusion through the first half of the channel to the "entry" water cluster near His37 may be hampered by significant kinetic barriers due to electrostatic repulsion. However, once in the entry cluster, a proton can move to one of the acceptor His37 in a nearly barrierless fashion, as evidenced both by MD simulations and a scan of the potential energy surface (PES). Water molecules of the entry cluster, although confined in the M2 pore and restricted in their motions, can conduct protons with a rate very similar to that of bulk water.

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.

On the role of water in peroxidase catalysis: a theoretical investigation of HRP compound I formation

We have investigated the dynamics of water molecules in the distal pocket of horseradish peroxidase to elucidate the role that they may play in the formation of the principal active species of the enzymatic cycle (compound I, Por(o+)-Fe(IV)=O) upon reaction of the resting Fe(III) state with hydrogen peroxide. The equilibrium molecular dynamics simulations show that, in accord with experimental evidence, the active site access channel is hydrated with an average of two to three water molecules within 5 A from the bound hydrogen peroxide. Although the channel is always hydrated, the specific conformations in which a water molecule bridges H(2)O(2) and the distal histidine, which were found (Derat; et al. J. Am. Chem. Soc. 2007, 129, 6346.) to display a low-energy barrier for the initial acid-base step of the reaction, occur with low probability but are relevant within the time scale of catalysis. Metadynamics simulations, which were used to reconstruct the free-energy landscape of water motion in the access channel, revealed that preferred interaction sites within the channel are separated by small energy barriers (<1.5 kcal/mol). Most importantly, water-bridged conformations lie on a shoulder just 1 kcal/mol above one local minimum and thus are easily accessible. Such an energy landscape appears as a requisite for the effectiveness of compound I formation, whereby the H-bonding pattern involving reactants and catalytic residues (including the intervening water molecule) has to rearrange to deliver the proton to the distal OH moiety of the hydrogen peroxide and thereby lead to heterolytic O-O cleavage. Our study provides an example of a system for which the "reactive configurations" (i.e., structures characterized by a low barrier for the chemical transformation) correspond to a minor population of the system and show how equilibrium molecular dynamics and free-energy calculations may conveniently be used to ascertain that such reactive conformations are indeed accessible to the system. Once again, the MD and QM/MM combination shows that a single water molecule acts as a biocatalyst in the cycle of HRP.

Functional studies and modeling of pore-lining residue mutants of the influenza a virus M2 ion channel

The A/M2 protein of influenza A virus forms a tetrameric proton-selective pH-gated ion channel. The H(37)xxxW(41) motif located in the channel pore is responsible for its gating and proton selectivity. Channel activation most likely involves protonation of the H37 residues, while the conductive state of the channel is characterized by two or three charged His residues in a tetrad. A/M2 channel activity is inhibited by the antiviral drug amantadine. Although a large number of functional amantadine-resistant mutants of A/M2 have been observed in vitro, only a few are observed in highly transmissible viruses in the presence or absence of amantadine. We therefore examined 49 point mutants of the pore-lining residues, representing both natural and nonnatural variants. Their ion selectivity, amantadine sensitivity, specific activity, and pH-dependent conductance were measured in Xenopus oocytes. These measurements showed how variations in the sequence lead to variations in the proton conduction. The results are consistent with a multistep mechanism that allows the protein to fine-tune its pH-rate profile over a wide range of proton concentrations, hypothesized to arise from different protonation states of the H37 tetrad. Mutations that give native-like conductance at low pH as well as minimal leakage current at pH 7.0 were surprisingly rare. Moreover, the results are consistent with a location of the amantadine-binding site inside the channel pore. These findings have helped to define the set of functionally fit mutants that should be targeted when considering the design of novel drugs that inhibit amantadine-resistant strains of influenza A virus.

Exploring Multidimensional Free Energy Landscapes Using Time-Dependent Biases on Collective Variables

A new implementation of the adaptive biasing force (ABF) method is described. This implementation supports a wide range of collective variables and can be applied to the computation of multidimensional energy profiles. It is provided to the community as part of a code that implements several analogous methods, including metadynamics. ABF and metadynamics have not previously been tested side by side on identical systems. Here, numerical tests are carried out on processes including conformational changes in model peptides and translocation of a halide ion across a lipid membrane through a peptide nanotube. On the basis of these examples, we discuss similarities and differences between the ABF and metadynamics schemes. Both approaches provide enhanced sampling and free energy profiles in quantitative agreement with each other in different applications. The method of choice depends on the dimension of the reaction coordinate space, the height of the barriers, and the relaxation times of degrees of freedom in the orthogonal space, which are not explicitly described by the chosen collective variables.

Using metadynamics to understand the mechanism of calmodulin/target recognition at atomic detail

The ability of calcium-bound calmodulin (CaM) to recognize most of its target peptides is caused by its binding to two hydrophobic residues ('anchors'). In most of the CaM complexes, the anchors pack against the hydrophobic pockets of the CaM domains and are surrounded by fully conserved Met side chains. Here, by using metadynamics simulations, we investigate quantitatively the energetics of the final step of this process using the M13 peptide, which has a high affinity and spans the sequence of the skeletal myosin light chain kinase, an important natural CaM target. We established the accuracy of our calculations by a comparison between calculated and NMR-derived structural and dynamical properties. Our calculations provide novel insights into the mechanism of protein/peptide recognition: we show that the process is associated with a free energy gain similar to that experimentally measured for the CaM complex with the homologous smooth muscle MLCK peptide (Ehrhardt et al., 1995, Biochemistry 34, 2731). We suggest that binding is dominated by the entropic effect, in agreement with previous proposals. Furthermore, we explain the role of conserved methionines by showing that the large flexibility of these side chains is a key feature of the binding mechanism. Finally, we provide a rationale for the experimental observation that in all CaM complexes the C-terminal domain seems to be hierarchically more important in establishing the interaction.

Unwinding the helical linker of calcium-loaded calmodulin: A molecular dynamics study

The fold of calmodulin (CaM) consists of two globular domains connected by a helical segment (the linker), whose conformational properties play a crucial role for the protein's molecular recognition processes. Here we investigate the structural properties of the linker by performing a 11.5 ns molecular dynamics (MD) simulation of calcium-loaded human CaM in aqueous solution. The calculations are based on the AMBER force field. The calculated S2 order parameters are in good accord with NMR data: The structure of the linker in our simulations is much more flexible than that emerging from the Homo sapiens X-ray structure, consistently with the helix unwinding observed experimentally in solution. This process occurs spontaneously in a nanosecond timescale, as observed also in a very recent simulation based on the GROMOS force field. A detailed description of the mechanism that determines the linker unwinding is provided, in which electrostatic contacts between the two globular domains play a critical role. The orientation of the domains emerging from our MD calculations is consistent both with former X-ray scattering data and a recent NMR work. Based on our findings, a rationale for the experimentally measured entropy cost associated to binding to the protein's cellular partners is also given.