Molecular dynamics simulations of membrane proteins

Membrane proteins control the traffic across cell membranes and thereby play an essential role in cell function from transport of various solutes to immune response via molecular recognition. Because it is very difficult to determine the structures of membrane proteins experimentally, computational methods have been increasingly used to study their structure and function. Here we focus on two classes of membrane proteins—ion channels and transporters—which are responsible for the generation of action potentials in nerves, muscles, and other excitable cells. We describe how computational methods have been used to construct models for these proteins and to study the transport mechanism. The main computational tool is the molecular dynamics (MD) simulation, which can be used for everything from refinement of protein structures to free energy calculations of transport processes. We illustrate with specific examples from gramicidin and potassium channels and aspartate transporters how the function of these membrane proteins can be investigated using MD simulations.     

Test of molecular dynamics force fields in gramicidin A

The force fields commonly used in molecular dynamics simulations of proteins are optimized under bulk conditions. Whether the same force fields can be used in simulations of membrane proteins is not well established, although they are increasingly being used for such purposes. Here we consider ion permeation in the gramicidin A channel as a test of the AMBER force field in a membrane environment. The potentials of mean force for potassium ions are calculated along the channel axis and compared with the one deduced from the experimental conductance data. The calculated result indicates a rather large central barrier similar to those obtained from other force fields, which are incompatible with the conductance data. We suggest that lack of polarizability is the most likely cause of this problem, and, therefore, urge development of polarizable force fields for simulations of membrane proteins.     

Modelling of β-d-glucopyranose ring distortion in different force fields: a metadynamics study

Modelling of carbohydrate conformations is a challenging task for force field developers. Three carbohydrate force fields, namely GLYCAM06, GROMOS 45a4 and OPLS were evaluated. Free energies of different ring conformations of β-d-glucopyranose were calculated using metadynamics in vacuum as well as in explicitly modelled water. All three force fields model the ⁴C₁ conformation as the most stable by at least 6 kJ/mol, as compared to other conformations. Interconversion from the ⁴C₁ to any other conformation is associated with a barrier of no lower than 26 kJ/mol. The free energy surface calculated in the GLYCAM06 force field is in remarkably good agreement with the recent Car-Parrinello metadynamics study. The effect of a water environment is relatively low and analogous in all tested force fields. Namely, the presence of water stabilizes the upper-left (^3,OB) versus bottom-right (B_3,O) area of Stoddard’s plot, relative to the situation in vacuum. Comparison of free and potential surfaces is also provided for vacuum calculations.     

Molecular-dynamics simulations of pyronine 6G and rhodamine 6G dimers in aqueous solution

We have carried out molecular-dynamics (MD) simulations on dimers of the positively charged laser dyes pyronine 6G (P6G) and rhodamine 6G (R6G) in aqueous solution, generating trajectories of 2.5 ns for various computational protocols. We discuss how the choice of atomic partial charges and the length of the trajectories affect the predicted structures of the dimers and compare our results to those of earlier MD-simulations, which were restricted to only 0.7 ns. Our results confirm that monomers of P6G easily undergo relative rotations within the dimer, but we found new conformations of the R6G dimer at longer simulation times. In addition, we analyzed in detail the energy change during the formation of dimers. With suitable corrections, the electrostatic energy from an Ewald treatment agrees with the results from an approach relying on a residue-based cutoff. For P6G, we show that the strong solvent-mediated electrostatic attraction between the monomers is counteracted by an almost equally large solvent-induced entropy contribution to yield a small driving force to dimer formation, in very good agreement with the free-energy change from a thermodynamic-integration procedure. Thus, earlier rationalizations of the dimer formation, based only on energy arguments, yield a qualitatively wrong picture.Electronic Supplementary Material Supplementary material is available for this article at and is accessible for authorized users.     

Molecular dynamics - potential of mean force calculations as a tool for understanding ion permeation and selectivity in narrow channels

Ion channels catalyze the permeation of charged molecules across cell membranes and are essential for many vital physiological functions, including nerve and muscle activity. To understand better the mechanisms underlying ion conduction and valence selectivity of narrow ion channels, we have employed free energy techniques to calculate the potential of mean force (PMF) for ion movement through the prototypical gramicidin A channel. Employing modern all-atom molecular dynamics (MD) force fields with umbrella sampling methods that incorporate one hundred 1-2 ns trajectories, we find that it is possible to achieve semi-quantitative agreement with experimental binding and conductance measurements. We also examine the sensitivity of the MD-PMF results to the choice of MD force field and compare PMFs for potassium, calcium and chloride ions to explore the basis for the valence selectivity of this narrow and uncharged ion channel. A large central barrier is observed for both anions and divalent ions, consistent with lack of experimental conductance. Neither anion or divalent cation is seen to be stabilized inside the channel relative to the bulk electrolyte and each leads to large disruptions to the protein and membrane structure when held deep inside the channel. Weak binding of calcium ions outside the channel corresponds to a free energy well that is too shallow to demonstrate channel blocking. Our findings emphasize the success of the MD-PMF approach and the sensitivity of ion energetics to the choice of biomolecular force field.     

KcsA closed and open: modelling and simulation studies

Bacterial homologues of mammalian potassium channels provide structures of two states of a gated K channel. Thus, the crystal structure of KcsA represents a closed state whilst that of MthK represents an open state. Using homology modelling and molecular dynamics simulations we have built a model of the transmembrane domain of KcsA in an open state and have compared its conformational stability with that of the same domain of KcsA in a closed state. Approximate Born energy calculations of monovalent cations within the two KcsA channel states suggest that the intracellular hydrophobic gate in the closed state provides a barrier of height ~5 kT to ion permeation, whilst in the open state the barrier is absent. Simulations (10 ns duration) in an octane slab (a simple membrane mimetic) suggest that closed- and open-state models are of comparable conformational stability, both exhibiting conformational drifts of ~3.3 Å C RMSD relative to the respective starting models. Substantial conformational fluctuations are observed in the intracellular gate region during both simulations (closed state and open state). In the simulation of open-state KcsA, rapid (<5 ns) exit of all three K⁺ ions occurs through the intracellular mouth of the channel. Helix kink and swivel motion is observed at the molecular hinge formed by residue G99 of the M2 helix. This motion is more substantial for the open- than for the closed-state model of the channel.     

Exploring blocker binding to a homology model of the open hERG K+ channel using docking and molecular dynamics methods

Binding of blockers to the human voltage-gated hERG potassium channel is studied using a combination of homology modelling, automated docking calculations and molecular dynamics simulations, where binding affinities are evaluated using the linear interaction energy method. A homology model was constructed based on the available crystal structure of the bacterial KvAP channel and the affinities of a series of sertindole analogues predicted using this model. The calculations reproduce the relative binding affinities of these compounds very well and indicate that both polar interactions near the intracellular opening of the selectivity filter as well as hydrophobic complementarity in the region around F656 are important for blocker binding. These results are consistent with recent alanine scanning mutation experiments on the blocking of the hERG channel by other compounds.     

Selection and flexible optimization of binding modes from conformation ensembles

This paper describes a new semi-flexible docking approach named Fado (flexible alignment and docking), which incorporates flexibility by using an ensemble of precomputed ligand conformers. A primary ligand is defined as a linear combination over all input conformers. An optimization with regard to the linear coefficients makes the ligand flexible. Initially, a point matching problem utilizing the Merck Molecular Force Field (MMFF) is modeled in order to compute the correct orientation of the ligand with respect to the target. The problem is then solved through a local optimization approach (RPROP). This is done for 20 randomized ligand orientations, yielding 20 binding modes per complex. Evaluating these modes illustrates that our method is able to reproduce the binding modes of molecules within a few minutes of CPU time. A representative dataset of diverse protein-ligand complexes could be reproduced with 78% accuracy below 2A RMSD distance to the reference crystal structure. Fado is available upon request to the authors (see also http://www.zib.eu/Numerik/projects/docking/projectlong.en.html).     

A systematic molecular dynamics simulation study of temperature dependent bilayer structural properties

Although lipid force fields (FFs) used in molecular dynamics (MD) simulations have proved to be accurate, there has not been a systematic study on their accuracy over a range of temperatures. Motivated by the X-ray and neutron scattering measurements of common phosphatidylcholine (PC) bilayers (Ku?erka et al. BBA. 1808: 2761, 2011), the CHARMM36 (C36) FF accuracy is tested in this work with MD simulations of six common PC lipid bilayers over a wide range of temperatures. The calculated scattering form factors and deuterium order parameters from the C36 MD simulations agree well with the X-ray, neutron, and NMR experimental data. There is excellent agreement between MD simulations and experimental estimates for the surface area per lipid, bilayer thickness (D_B), hydrophobic thickness (D_C), and lipid volume (V_L). The only minor discrepancy between simulation and experiment is a measure of (D_B − D_HH) / 2 where D_HH is the distance between the maxima in the electron density profile along the bilayer normal. Additional MD simulations with pure water and heptane over a range of temperatures provide explanations of possible reasons causing the minor deviation. Overall, the C36 FF is accurate for use with liquid crystalline PC bilayers of varying chain types and over biologically relevant temperatures.     

Reevaluation of fenpropimorph as a σ receptor ligand: Structure-affinity relationship studies at human σ1 receptors

Fenpropimorph (1) is considered a “super high-affinity” σ₁ receptor ligand (K_i = 0.005 nM for guinea pig σ₁ receptors). Here, we examine the binding of 1 and several of its deconstructed analogs at human σ₁ (hσ₁) receptors. We monitored their subtype selectivity by determining the binding affinity at σ₂ receptors. In addition, we validated an existing pharmacophore model at the molecular level by conducting 3D molecular modeling studies, using the crystal structure of hσ₁ receptors, and Hydrophatic INTeractions (HINT) analysis. Our structure affinity relationship studies showed that 1 binds with lower affinity at hσ₁ receptors (K_i = 17.3 nM) compared to guinea pig; moreover, we found that none of the fenpropimorph methyl groups is important for its binding at hσ₁ receptors, nor is stereochemistry. For example, removal of all methyl groups as seen in 4 resulted in an almost 5-fold higher affinity at hσ₁ receptors compared to 1 and 350-fold selectivity versus σ₂ receptors. In addition, although the O atom of the morpholine ring does not contribute to affinity at hσ₁ receptors (and might even detract from it), it plays role in subtype (σ₁ versus σ₂ receptor) selectivity.     

Influence of protein flexibility on the electrostatic energy landscape in gramicidin A

We describe an electrostatic model of the gramicidin A channel that allows protein atoms to move in response to the presence of a permeating ion. To do this, molecular dynamics simulations are carried out with a permeating ion at various positions within the channel. Then an ensemble of atomic coordinates taken from the simulations are used to construct energy profiles using macroscopic electrostatic calculations. The energy profiles constructed are compared to experimentally-determined conductance data by inserting them into Brownian dynamics simulations. We find that the energy landscape seen by a permeating ion changes significantly when we allow the protein atoms to move rather than using a rigid protein structure. However, the model developed cannot satisfactorily reproduce all of the experimental data. Thus, even when protein atoms are allowed to move, the dielectric model used in our electrostatic calculations breaks down when modeling the gramicidin channel.     

Continuous metadynamics in essential coordinates as a tool for free energy modelling of conformational changes

Modelling of conformational changes in biopolymers is one of the greatest challenges of molecular biophysics. Metadynamics is a recently introduced free energy modelling technique that enhances sampling of configurational (e.g. conformational) space within a molecular dynamics simulation. This enhancement is achieved by the addition of a history-dependent bias potential, which drives the system from previously visited regions. Discontinuous metadynamics in the space of essential dynamics eigenvectors (collective motions) has been proposed and tested in conformational change modelling. Here, we present an implementation of two continuous formulations of metadynamics in the essential subspace. The method was performed in a modified version of the molecular dynamics package GROMACS. These implementations were tested on conformational changes in cyclohexane, alanine dipeptide (terminally blocked alanine, Ace-Ala-Nme) and SH3 domain. The results illustrate that metadynamics in the space of essential coordinates can accurately model free energy surfaces associated with conformational changes. Figure The conformational free energy surface of cyclohexane in the space of the two most intensive collective motions.

---

     

Fluctuation driven transport and models of molecular motors and pumps

Non-equilibrium fluctuations can drive vectorial transport along an anisotropic structure in an isothermal medium by biasing the effect of thermal noise (k _B T). Mechanisms based on this principle are often called Brownian ratchets and have been invoked as a possible explanation for the operation of biomolecular motors and pumps. We discuss the thermodynamics and kinetics for the operation of microscopic ratchet motors under conditions relevant to biology, showing how energy provided by external fluctuations or a non-equilibrium chemical reaction can cause unidirectional motion or uphill pumping of a substance. Our analysis suggests that molecular pumps such as Na,K-ATPase and molecular motors such as kinesin and myosin may share a common underlying mechanism. Received: 18 February 1998 / Revised version: 5 May 1998 / Accepted: 14 May 1998     

Free energy analysis of conductivity and charge selectivity of M2GlyR-derived synthetic channels

Significant progresses have been made in the design, synthesis, modeling and in vitro testing of channel-forming peptides derived from the second transmembrane domain of the α-subunit of the glycine receptor (GlyR). The latest designs, including p22 (KKKKP ARVGL GITTV LTMTT QS), are highly soluble in water with minimal aggregation propensity and insert efficiently into cell membranes to form highly conductive ion channels. The last obstacle to a potential lead sequence for channel replacement treatment of CF patients is achieving adequate chloride selectivity. We have performed free energy simulation to analyze the conductance and charge selectivity of M2GlyR-derived synthetic channels. The results reveal that the pentameric p22 pore is non-selective. Moderate barriers for permeation of both K⁺ and Cl⁻ are dominated by the desolvation cost. Despite previous evidence suggesting a potential role of threonine side chains in anion selectivity, the hydroxyl group is not a good surrogate of water for coordinating these ions. We have also tested initial ideas of introducing additional rings of positive changes to various positions along the pore to increase anion selectivity. The results support the feasibility of achieving anion selectivity by modifying the electrostatic properties of the pore, but at the same time suggest that the peptide assembly and pore topology may also be dramatically modified, which could abolish the effects of modified electrostatics on anion selectivity. This was confirmed by subsequent two-electrode voltage clamp measurements showing that none of the tested mono-, di- and tri-Dap substituted sequences was selective. The current study thus highlights the importance of controlling channel topology besides modifying pore electrostatics for achieving anion selectivity. Several strategies are now being explored in our continued efforts to design an anion selective peptide channel with suitable biophysical, physiological and pharmacological properties as a potential treatment modality for channel replacement therapy. This article is part of a Special Issue entitled: Interfacially Active Peptides and Proteins. Guest Editors: William C. Wimley and Kalina Hristova.     

生物大分子的自由能计算方法

在本文中,我们总结了计算生物物理中广泛使用的经典的自由能计算方法,主要可以归纳为三大类:第一类是基于过渡态理论,利用路径系综得到自由能的方法;第二类方法采用一些经验势来补偿自由能表面达到均匀采样的目标;最后一种是将整个转换径分段计算,然后拼接完成得到完整自由能的方法.通过详细介绍自由能计算的重要性和实施细节,希望能够让...