Insights into lipid-protein interactions from computer simulations.

Lipid-protein interactions play an important direct role in the function of many membrane proteins. We argue they are key players in membrane structure, modulate membrane proteins in more subtle ways than direct binding, and are important for understanding the mechanism of classes of hydrophobic drugs. By directly comparing membrane proteins from different families in the same, complex lipid mixture, we found a unique lipid environment for every protein. Extending this work, we identified both differences and similarities in the lipid environment of GPCRs, dependent on which family they belong to and in some cases their conformational state, with particular emphasis on the distribution of cholesterol. More recently, we have been studying modes of coupling between protein conformation and local membrane properties using model proteins. In more applied approaches, we have used similar methods to investigate specific hypotheses on interactions of lipid and lipid-like molecules with ion channels. We conclude this perspective with some considerations for future work, including a new more sophisticated coarse-grained force field (Martini 3), an interactive visual exploration framework, and opportunities to improve sampling.

Structural Properties of Inverted Hexagonal Phase: A Hybrid Computational and Experimental Approach.

Inverted/reverse hexagonal (HII) phases are of special interest in several fields of research, including nanomedicine. We used molecular dynamics (MD) simulation to study HII systems composed of dioleoylphosphatidylethanolamine (DOPE) and palmitoyloleoylphosphatidylethanolamine (POPE) at several hydration levels and temperatures. The effect of the hydration level on several HII structural parameters, including deuterium order parameters, was investigated. We further used MD simulations to estimate the maximum hydrations of DOPE and POPE HII lattices at several given temperatures. Finally, the effect of acyl chain unsaturation degree on the HII structure was studied via comparing the DOPE with POPE HII systems. In addition to MD simulations, we used deuterium nuclear magnetic resonance (2H NMR) and small-angle X-ray scattering (SAXS) experiments to measure the DOPE acyl chain order parameters, lattice plane distances, and the water core radius in HII phase DOPE samples at several temperatures in the presence of excess water. Structural parameters calculated from MD simulations are in excellent agreement with the experimental data. Dehydration decreases the radius of the water core. An increase in hydration level slightly increased the deuterium order parameter of lipids acyl chains, whereas an increase in temperature decreased it. Lipid cylinders undulated along the cylinder axis as a function of hydration level. The maximum hydration levels of PE HII phases at different temperatures were successfully predicted by MD simulations based on a single experimental measurement for the lattice plane distance in the presence of excess water. An increase in temperature decreases the maximum hydration and consequently the radius of the water core and lattice plane distances. Finally, DOPE formed HII structures with a higher curvature compared to POPE, as expected. We propose a general protocol for constructing computational HII systems that correspond to the experimental systems. This protocol could be used to study HII systems composed of molecules other than the PE systems used here and to improve and validate force field parameters by using the target data in the HII phase.

Molecular Ordering in Lipid Monolayers: An Atomistic Simulation.

We report on atomistic simulations of DPPC lipid monolayers using the CHARMM36 lipid force field (and also the Slipid force field as a control case), combined with a four-point OPC water model. The entire two-phase region where domains of the "liquid-condensed" (LC) phase coexist with domains of the "liquid-expanded" (LE) phase has been explored. The simulations are long enough that the complete phase-transition stage, with two domains coexisting in the monolayer, is reached in all cases. Also, system sizes used are larger than those in previous works. As expected, domains of the minority phase are elongated, emphasizing the importance of anisotropic van der Waals and/or electrostatic dipolar interactions in the monolayer plane. The molecular structure is quantified in terms of distribution functions for the hydrocarbon chains and the PN dipoles. In contrast to previous work, where average distributions are calculated, distributions are here extracted for each of the coexisting phases by first identifying lipid molecules that belong to either LC or LE regions. In the case of the CHARMM36 force field, the three-dimensional distributions show that the average tilt angle of the chains with respect to the normal outward direction is (39.0 ± 0.1)° in the LC phase and (48.1 ± 0.5)° in the LC phase. In the case of the PN dipoles, the distributions indicate a tilt angle of (110.8 ± 0.5)° in the LC phase and (112.5 ± 0.5)° in the LE phase. These results are quantitatively different from those in previous works, which indicated a smaller normal component of the PN dipole. Also, the distributions of the monolayer-projected chains and PN dipoles have been calculated. Chain distributions peak along a particular direction in the LC domains, while they are uniform in the LE phase. Long-range ordering associated with the projected PN dipoles is absent in both phases. These results strongly suggest that LC domains do not exhibit dipolar ordering in the plane of the monolayer, the effect of these components being averaged out at short distances. Therefore, the only relevant component of the molecular dipoles, with regard to both intra- and long-range interdomain interactions, is normal to the monolayer. Also, the local orientation of chain projections is almost constant in LC domains and points in the direction along which domains are elongated, suggesting that the line tension driving the phase transition might be anisotropic with respect to the interfacial domain boundary.

Ionizable amino lipid interactions with POPC: implications for lipid nanoparticle function.

Lipid nanoparticles (LNPs) composed of ionizable cationic lipids are currently the leading systems for siRNA delivery in liver disease, with the major limitation of low siRNA release efficacy into the cytoplasm. Ionizable cationic lipids are known to be of critical importance in LNP structure and stability, siRNA entrapment, and endosomal disruption. However, their distribution inside the LNPs and their exact role in cytoplasmic delivery remain unclear. A recent study [Kulkarni et al., On the formation and morphology of lipid nanoparticles containing ionizable cationic lipids and siRNA, ACS Nano, 2018, 12(5), 4787-4795] on LNP-siRNA systems containing the ionizable lipid DLin-KC2-DMA (also known as KC2 with an apparent pKa of ca. 6.7) suggested that neutral KC2 segregates from other components and forms an amorphous oil droplet in the core of LNPs. In this paper, we present evidence supporting the model proposed by Kulkarni et al. We studied KC2 segregation in the presence of POPC using molecular dynamics simulation, deuterium NMR, SAXS, and cryo-TEM experiments, and found that neutral KC2 has a high tendency to separate from POPC dispersions. KC2 confinement, upon raising the pH during the formulation process, could result in rearrangement of the internal structure of LNPs. As interactions between cationic KC2 and anionic endosomal lipids are thought to be a key factor in cargo release, KC2 confinement inside the LNP may be responsible for the observed low release efficacy.

Computational and experimental approaches for investigating nanoparticle-based drug delivery systems.

Most therapeutic agents suffer from poor solubility, rapid clearance from the blood stream, a lack of targeting, and often poor translocation ability across cell membranes. Drug/gene delivery systems (DDSs) are capable of overcoming some of these barriers to enhance delivery of drugs to their right place of action, e.g. inside cancer cells. In this review, we focus on nanoparticles as DDSs. Complementary experimental and computational studies have enhanced our understanding of the mechanism of action of nanocarriers and their underlying interactions with drugs, biomembranes and other biological molecules. We review key biophysical aspects of DDSs and discuss how computer modeling can assist in rational design of DDSs with improved and optimized properties. We summarize commonly used experimental techniques for the study of DDSs. Then we review computational studies for several major categories of nanocarriers, including dendrimers and dendrons, polymer-, peptide-, nucleic acid-, lipid-, and carbon-based DDSs, and gold nanoparticles. This article is part of a Special Issue entitled: Membrane Proteins edited by J.C. Gumbart and Sergei Noskov.

Interactions of key charged residues contributing to selective block of neuronal sodium channels by µ-conotoxin KIIIA.

Voltage-gated sodium channels are important in initiating and propagating nerve impulses in various tissues, including cardiac muscle, skeletal muscle, the brain, and the peripheral nerves. Hyperexcitability of these channels leads to such disorders as cardiac arrhythmias (Na(v)1.5), myotonias (Na(v)1.4), epilepsies (Na(v)1.2), and pain (Na(v)1.7). Thus, there is strong motivation to identify isoform-specific blockers and the molecular determinants underlying their selectivity among these channels. µ-Conotoxin KIIIA blocks rNa(v)1.2 (IC(50), 5 nM), rNa(v)1.4 (37 nM), and hNa(v)1.7 (97 nM), expressed in mammalian cells, with high affinity and a maximal block at saturating concentrations of 90 to 95%. Mutations of charged residues on both the toxin and channel modulate the maximal block and/or affinity of KIIIA. Two toxin substitutions, K7A and R10A, modulate the maximal block (52-70%). KIIIA-H12A and R14A were the only derivatives tested that altered Na(v) isoform specificity. KIIIA-R14A showed the highest affinity for Na(v)1.7, a channel involved in pain signaling. Wild-type KIIIA has a 2-fold higher affinity for Na(v)1.4 than for Na(v)1.7, which can be attributed to a missing outer vestibule charge in domain III of Na(v)1.7. Reciprocal mutations Na(v)1.4 D1241I and Na(v)1.7 I1410D remove the affinity differences between these two channels for wild-type KIIIA without affecting their affinities for KIIIA-R14A. KIIIA is the first µ-conotoxin to show enhanced activity as pH is lowered, apparently resulting from titration of the free N terminus. Removal of this free amino group reduced the pH sensitivity by 10-fold. Recognition of these molecular determinants of KIIIA block may facilitate further development of subtype-specific, sodium channel blockers to treat hyperexcitability disorders.

Orientation of µ-conotoxin PIIIA in a sodium channel vestibule, based on voltage dependence of its binding.

Mutant cycle analysis has been used in previous studies to constrain possible docking orientations for various toxins. As an independent test of the bound orientation of µ-conotoxin PIIIA, a selectively targeted sodium channel pore blocker, we determined the contributions to binding voltage dependence of specific residues on the surface of the toxin. A change in the "apparent valence" (zδ) of the block, which is associated with a change of a specific toxin charge, reflects a change in the charge movement within the transmembrane electric field as the toxin binds. Toxin derivatives with charge-conserving mutations (R12K, R14K, and K17R) showed zδ values similar to those of wild type (0.61 ± 0.01, mean ± S.E.M.). Charge-changing mutations produced a range of responses. Neutralizing substitutions for Arg14 and Lys17 showed the largest reductions in zδ values, to 0.18 ± 0.06 and 0.20 ± 0.06, respectively, whereas unit charge-changing substitutions for Arg12, Ser13, and Arg20 gave intermediate values (0.24 ± 0.07, 0.33 ± 0.04, and 0.32 ± 0.05), which suggests that each of these residues contributes to the dependence of binding on the transmembrane voltage. Two mutations, R2A and G6K, yielded no significant change in zδ. These observations suggest that the toxin binds with Arg2 and Gly6 facing the extracellular solution, and Arg14 and Lys17 positioned most deeply in the pore. In this study, we used molecular dynamics to simulate toxin docking and performed Poisson-Boltzmann calculations to estimate the changes in local electrostatic potential when individual charges were substituted on the toxin's surface. Consideration of two limiting possibilities suggests that most of the charge movement associated with toxin binding reflects sodium redistribution within the narrow part of the pore.

Structure and dynamics of the antifungal molecules Syringotoxin-B and Syringopeptin-25A from molecular dynamics simulation.

The phytopathogen Pseudomonas syringae pv. syringae produces toxic cyclic lipodepsipeptides (CLPs): nona-peptides and syringopeptins. All CLPs inhibit the growth of many fungal species, including human pathogens, although different fungi display different degrees of sensitivity. The best studied CLPs are Syringomycin-E (SR-E), Syringotoxin-B (ST-B) and Syringopeptin-25A (SP-25A). Their biological activity is affected by membrane composition and their structural differences. We previously (Matyus et al. in Eur Biophys J 35:459-467, 2006) reported the molecular features and structural preferences of SR-E in water and octane environments. Here we investigate in atomic detail the molecular features of the two other main CLP components, ST-B and SP-25A, in water and octane by 200 ns molecular dynamics simulations (MD), using distance restraints derived from NMR NOE data (Ballio et al. in Eur J Biochem 234:747-758, 1995). We have obtained three-dimensional models of ST-B and SP-25A CLPs in different environments. These models can now be used as a basis to investigate the interactions of ST-B and SP-25A with lipid membranes an important further step towards a better understanding of the antifungal and antibacterial activity of these peptides.

Structural investigation of syringomycin-E using molecular dynamics simulation and NMR.

Syringomycin-E (SR-E) is a cyclic lipodepsinonapeptide produced by certain strains of the bacterium Pseudomonas syringae pv. syringae. It shows inhibitory effects against many fungal species, including human pathogens. Its primary biological target is the plasma membrane, where it forms channels comprised of at least six SR-E molecules. The high-resolution structure of SR-E and the structure of the channels are currently not known. In this paper, we investigate in atomic detail the molecular features of SR-E in water by NMR and in water and octane by molecular dynamics simulation (MD). We built a model of the peptide and examined its structure in water and octane in 200 ns MD simulations both with and without distance restraints derived from NMR NOE data. The resulting trajectories show good agreement with the measured NOEs and circular dichroism data from the literature and provide atomistic models of SR-E that are an important step toward a better understanding of the antifungal and antibacterial activity of this peptide.

Orientation and dynamics of benzyl alcohol and benzyl alkyl ethers dissolved in nematic lyotropic liquid crystals. 2H NMR and molecular dynamics simulations.

Most drugs have to cross cell membranes to reach their final target. A better understanding of the distribution, interactions, and dynamics of biologically active molecules in model bilayers is of fundamental importance in understanding drug functioning and design. 2H NMR quadrupole splittings (delta nu(Q)) and longitudinal relaxation times (T1) from the aromatic ring of benzyl alcohol-d5 (C0), a commonly used anesthetic, and a series of linear alkyl benzyl-d5 ethers with chain lengths from 1 to 12 carbon atoms (C1-C12), were measured. The molecules were dissolved in a nematic discotic lyotropic liquid crystal solution made of tetradecyltrimethylammonium chloride (TTAC)/decanol (DeOH)/NaCl/H2O. Values of delta nu(Q) and T1 from 1,1-dideuteriodecanol (15% enriched) and DHO (H2O with 0.2% D2O) were also measured. Delta nu(Q) of DeOH and DHO remained constant throughout the series. The value of delta nu(Q) of the para position of the ring (delta nu(p)) in C1 is 30% smaller than the delta nu(p) of C0. This is attributed to the existence of an H-bond between the alcohol hydroxyl proton and the solvent, which influences the average orientation of the ring. The relaxation data show that T1o,m is always longer than T1p and both decrease with the increase in alkyl chain length. Molecular dynamics simulations of the experimentally studied systems were performed. The aggregate was represented as a bilayer. The distribution, average orientation, and order parameters of the aromatic ring of the guest molecules in the bilayer were examined. Rotational correlation functions of all the C-D bonds and the OH bond from H2O were evaluated, allowing an estimate of the correlation times and T1. According to these results all spins relax in extreme narrowing conditions, except DeOH. Experimental and calculated T1 values differ at most by a factor of 3. However, the order of magnitude and the observed trends are well reproduced by the calculations. The aromatic ring of C0 possesses a unique average orientation in the bilayer. For the ether series, the orientation is modified and the C2 symmetry axis of the aromatic ring is exchanging between two orientations averaging the quadrupole splittings from the ortho and meta positions. The simulation supports the existence of an H-bond between C0 and the solvent not found in the ethers, which should be responsible for the observed differences.

Structure and aggregation number of a lyotropic liquid crystal: a fluorescence quenching and molecular dynamics study.

The structure and aggregation number of a discotic lyotropic liquid crystal, prepared from tetradecyltrimethylammonium chloride (TDTMACl)/decanol (DeOH)/NaCl/H2O, have been examined using fluorescence quenching of pyrene by hexadecylpyridinium chloride and molecular dynamics (MD). The fluorescence method gives an aggregation number of 258 +/- 25 units (DeOH + TDTMACl). From the MD simulation, a lower limit for the aggregate dimension of 130 units of DeOH + TDTMACl is predicted. A stable oblate aggregate of 240 units was studied in detail. A strong polarization between the ammonium headgroups and chloride ions is observed from the calculated trajectory. DeOH headgroups are located, on average, 0.3 nm more to the interior of the aggregate than the TDTMACl headgroup and contribute to widening the interface by forming H-bonds with water. The radial distribution function of the ammonium headgroup shows that there are 16 water molecules in the first solvation sphere. The diagonal elements of the order parameter tensor of the tail atoms of both surfactants indicate that the interior of the micelle preserves about the same degree of order as at the interface, up to the last three atoms of the aliphatic chain, where the order starts to decrease.

Pores formed by the nicotinic receptor m2delta Peptide: a molecular dynamics simulation study.

The M2delta peptide self-assembles to form a pentameric bundle of transmembrane alpha-helices that is a model of the pore-lining region of the nicotinic acetylcholine receptor. Long (>15 ns) molecular dynamics simulations of a model of the M2delta(5) bundle in a POPC bilayer have been used to explore the conformational dynamics of the channel assembly. On the timescale of the simulation, the bundle remains relatively stable, with the polar pore-lining side chains remaining exposed to the lumen of the channel. Fluctuations at the helix termini, and in the helix curvature, result in closing/opening transitions at both mouths of the channel, on a timescale of approximately 10 ns. On average, water within the pore lumen diffuses approximately 4x more slowly than water outside the channel. Examination of pore water trajectories reveals both single-file and path-crossing regimes to occur at different times within the simulation.

Molecular dynamics simulation of a lipid diamond cubic phase.

This paper presents the first atomistic simulation of a cubic membrane phase. Using the molecular dynamics simulation technique both the global and the local organization of glycerolmonoolein molecules inside the diamond cubic phase are studied. Multinanosecond simulations reveal that the center of the cubic bilayer remains close to the infinite periodic minimal surface that describes the diamond geometry. We further show that the equilibrium structure of the surfactant molecules inside the cubic phase is very similar to their structure inside a simulated lamellar bilayer. The small differences arise from the packing constraints of the surfactants within the cubic phase which has an area per surfactant that increases toward the bilayer center.

Proline-induced hinges in transmembrane helices: possible roles in ion channel gating.

A number of ion channels contain transmembrane (TM) alpha-helices that contain proline-induced molecular hinges. These TM helices include the channel-forming peptide alamethicin (Alm), the S6 helix from voltage-gated potassium (Kv) channels, and the D5 helix from voltage-gated chloride (CLC) channels. For both Alm and KvS6, experimental data implicate hinge-bending motions of the helix in an aspect of channel gating. We have compared the hinge-bending motions of these TM helices in bilayer-like environments by multi-nanosecond MD simulations in an attempt to describe motions of these helices that may underlie possible modes of channel gating. Alm is an alpha-helical channel-forming peptide, which contains a central kink associated with a Gly-x-x-Pro motif in its sequence. Simulations of Alm in a TM orientation for 10 ns in an octane slab indicate that the Gly-x-x-Pro motif acts as a molecular hinge. The S6 helix from Shaker Kv channels contains a Pro-Val-Pro motif. Modeling studies and recent experimental data suggest that the KvS6 helix may be kinked in the vicinity of this motif. Simulations (10 ns) of an isolated KvS6 helix in an octane slab and in a POPC bilayer reveal hinge-bending motions. A pattern-matching approach was used to search for possible hinge-bending motifs in the TM helices of other ion channel proteins. This uncovered a conserved Gly-x-Pro motif in TM helix D5 of CLC channels. MD simulations of a model of hCLC1-D5 spanning an octane slab suggest that this channel also contains a TM helix that undergoes hinge-bending motion. In conclusion, our simulations suggest a model in which hinge-bending motions of TM helices may play a functional role in the gating mechanisms of several different families of ion channels.

Voltage-dependent insertion of alamethicin at phospholipid/water and octane/water interfaces.

Understanding the binding and insertion of peptides in lipid bilayers is a prerequisite for understanding phenomena such as antimicrobial activity and membrane-protein folding. We describe molecular dynamics simulations of the antimicrobial peptide alamethicin in lipid/water and octane/water environments, taking into account an external electric field to mimic the membrane potential. At cis-positive potentials, alamethicin does not insert into a phospholipid bilayer in 10 ns of simulation, due to the slow dynamics of the peptide and lipids. However, in octane N-terminal insertion occurs at field strengths from 0.33 V/nm and higher, in simulations of up to 100 ns duration. Insertion of alamethicin occurs in two steps, corresponding to desolvation of the Gln7 side chain, and the backbone of Aib10 and Gly11. The proline induced helix kink angle does not change significantly during insertion. Polyalanine and alamethicin form stable helices both when inserted in octane and at the water/octane interface, where they partition in the same location. In water, both polyalanine and alamethicin partially unfold in multiple simulations. We present a detailed analysis of the insertion of alamethicin into the octane slab and the influence of the external field on the peptide structure. Our findings give new insight into the mechanism of channel formation by alamethicin and the structure and dynamics of membrane-associated helices.

Towards the design and computational characterization of a membrane protein.

The design of a transmembrane four-helix bundle is described. We start with an idealized four-helix bundle geometry, then use statistical information to build a plausible transmembrane bundle. Appropriate residues are chosen using database knowledge on the sequences of membrane helices and loops, then the packing of the bundle core is optimized, and favorable side chain rotamers from rotamer libraries are selected. Next, we use explicit physical knowledge from biomolecular simulation force fields and molecular dynamics simulations to test whether the designed structure is physically possible. These procedures test whether the designed protein will indeed be alpha-helical, well packed and stable over a time scale of several nanoseconds in a realistic lipid bilayer environment. We then test a modeling approach that does not include sophisticated database knowledge about proteins, but rather relies on applying our knowledge of the physics that governs protein motions. This independent validation of the design is based on simulated annealing and restrained molecular dynamics simulation in vacuo, comparable to procedures used to refine NMR and X-ray structures.

Structure and dynamics of the pore-lining helix of the nicotinic receptor: MD simulations in water, lipid bilayers, and transbilayer bundles.

Multiple nanosecond duration molecular dynamics simulations on the pore-lining M2 helix of the nicotinic acetylcholine receptor reveal how its structure and dynamics change as a function of environment. In water, the M2 helix partially unfolds to form a molecular hinge in the vicinity of a central Leu residue that has been implicated in the mechanism of ion channel gating. In a phospholipid bilayer, either as a single transmembrane helix, or as part of a pentameric helix bundle, the M2 helix shows less flexibility, but still exhibits a kink in the vicinity of the central Leu. The single M2 helix tilts relative to the bilayer normal by 12 degrees, in agreement with recent solid state NMR data (Opella et al., Nat Struct Biol 6:374-379, 1999). The pentameric helix bundle, a model for the pore domain of the nicotinic receptor and for channels formed by M2 peptides in a bilayer, is remarkably stable over a 2-ns MD simulation in a bilayer, provided one adjusts the pK(A)s of ionizable residues to their calculated values (when taking their environment into account) before starting the simulation. The resultant transbilayer pore shows fluctuations at either mouth which transiently close the channel. Proteins 2000;39:47-55.

Exploring models of the influenza A M2 channel: MD simulations in a phospholipid bilayer.

The M2 protein of influenza A virus forms homotetrameric helix bundles, which function as proton-selective channels. The native form of the protein is 97 residues long, although peptides representing the transmembrane section display ion channel activity, which (like the native channel) is blocked by the antiviral drug amantadine. As a small ion channel, M2 may provide useful insights into more complex channel systems. Models of tetrameric bundles of helices containing either 18 or 22 residues have been simulated while embedded in a fully hydrated 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphatidylcholine bilayer. Several different starting models have been used. These suggest that the simulation results, at least on a nanosecond time scale, are sensitive to the exact starting structure. Electrostatics calculations carried out on a ring of four ionizable aspartate residues at the N-terminal mouth of the channel suggest that at any one time, only one will be in a charged state. Helix bundle models were mostly stable over the duration of the simulation, and their helices remained tilted relative to the bilayer normal. The M2 helix bundles form closed channels that undergo breathing motions, alternating between a tetramer and a dimer-of-dimers structure. Under these conditions either the channel forms a pocket of trapped waters or it contains a column of waters broken predominantly at the C-terminal mouth of the pore. These waters exhibit restricted motion in the pore and are effectively "frozen" in a way similar to those seen in previous simulations of a proton channel formed by a four-helix bundle of a synthetic leucine-serine peptide (, Biophys. J. 77:2400-2410).