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.     

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.     

Water permeation through gramicidin A: desformylation and the double helix: a molecular dynamics study

Multinanosecond molecular dynamics simulations of gramicidin A embedded in a dimyristoylphosphatidylcholine bilayer show a remarkable structural stability for both experimentally determined conformations: the head-to-head helical dimer and the double helix. Water permeability was found to be much higher in the double helical conformation, which is explained by lower hydrogen bond-mediated enthalpic barriers at the channel entrance and its larger pore size. Free-energy perturbation calculations show that the double helical structure is stabilized by the positive charges at the N termini introduced by the desformylation, whereas the helical dimer is destabilized. Together with the recent experimental observation that desformyl gramicidin conducts water hundredfold better than gramicidin, this suggests that desformyl gramicidin A predominantly occurs in the double helical conformation.     

Improved model of hydrated calcium ion for molecular dynamics simulations using classical biomolecular force fields

Calcium ions (Ca²⁺) play key roles in various fundamental biological processes such as cell signaling and brain function. Molecular dynamics (MD) simulations have been used to study such interactions, however, the accuracy of the Ca²⁺ models provided by the standard MD force fields has not been rigorously tested. Here, we assess the performance of the Ca²⁺ models from the most popular classical force fields AMBER and CHARMM by computing the osmotic pressure of model compounds and the free energy of DNA–DNA interactions. In the simulations performed using the two standard models, Ca²⁺ ions are seen to form artificial clusters with chloride, acetate, and phosphate species; the osmotic pressure of CaAc₂ and CaCl₂ solutions is a small fraction of the experimental values for both force fields. Using the standard parameterization of Ca²⁺ ions in the simulations of Ca²⁺‐mediated DNA–DNA interactions leads to qualitatively wrong outcomes: both AMBER and CHARMM simulations suggest strong inter‐DNA attraction whereas, in experiment, DNA molecules repel one another. The artificial attraction of Ca²⁺ to DNA phosphate is strong enough to affect the direction of the electric field‐driven translocation of DNA through a solid‐state nanopore. To address these shortcomings of the standard Ca²⁺ model, we introduce a custom model of a hydrated Ca²⁺ ion and show that using our model brings the results of the above MD simulations in quantitative agreement with experiment. Our improved model of Ca²⁺ can be readily applied to MD simulations of various biomolecular systems, including nucleic acids, proteins and lipid bilayer membranes. © 2016 Wiley Periodicals, Inc. Biopolymers 105: 752–763, 2016.     

Ion transport in the gramicidin channel: molecular dynamics study of single and double occupancy.

The structural and thermodynamic factors responsible for the singly and doubly occupied saturation states of the gramicidin channel are investigated with molecular dynamics simulations and free energy perturbation methods. The relative free energy of binding of all of the five common cations Li+, Na+, K+, Rb+, and Cs+ is calculated in the singly and doubly occupied channel and in bulk water. The atomic system, which includes the gramicidin channel, a model membrane made of neutral Lennard-Jones particles and 190 explicit water molecules to form the bulk region, is similar to the one used in previous work to calculate the free energy profile of a Na+ ion along the axis of the channel. In all of the calculations, the ions are positioned in the main binding sites located near the entrances of the channel. The calculations reveal that the doubly occupied state is relatively more favorable for the larger ions. Thermodynamic decomposition is used to show that the origin of the trend observed in the calculations is due to the loss of favorable interactions between the ion and the single file water molecules inside the channel. Small ions are better solvated by the internal water molecules in the singly occupied state than in the doubly occupied state; bigger ions are solvated almost as well in both occupation states. Water-channel interactions play a role in the channel response. The observed trends are related to general thermodynamical properties of electrolyte solutions.     

Gramicidin A channel as a test ground for molecular dynamics force fields

We use the well-known structural and functional properties of the gramicidin A channel to test the appropriateness of force fields commonly used in molecular dynamics (MD) simulations of ion channels. For this purpose, the high-resolution structure of the gramicidin A dimer is embedded in a dimyristoylphosphatidylcholine bilayer, and the potential of mean force of a K(+) ion is calculated along the channel axis using the umbrella sampling method. Calculations are performed using two of the most common force fields in MD simulations: CHARMM and GROMACS. Both force fields lead to large central barriers for K(+) ion permeation, that are substantially higher than those deduced from the physiological data by inverse methods. In long MD simulations lasting over 60 ns, several ions are observed to enter the binding site but none of them crossed the channel despite the presence of a large driving field. The present results, taken together with many earlier studies, highlights the shortcomings of the standard force fields used in MD simulations of ion channels and calls for construction of more appropriate force fields for this purpose.     

Simulations of ion permeation through a potassium channel: molecular dynamics of KcsA in a phospholipid bilayer

Potassium channels enable K(+) ions to move passively across biological membranes. Multiple nanosecond-duration molecular dynamics simulations (total simulation time 5 ns) of a bacterial potassium channel (KcsA) embedded in a phospholipid bilayer reveal motions of ions, water, and protein. Comparison of simulations with and without K(+) ions indicate that the absence of ions destabilizes the structure of the selectivity filter. Within the selectivity filter, K(+) ions interact with the backbone (carbonyl) oxygens, and with the side-chain oxygen of T75. Concerted single-file motions of water molecules and K(+) ions within the selectivity filter of the channel occur on a 100-ps time scale. In a simulation with three K(+) ions (initially two in the filter and one in the cavity), the ion within the central cavity leaves the channel via its intracellular mouth after approximately 900 ps; within the cavity this ion interacts with the Ogamma atoms of two T107 side chains, revealing a favorable site within the otherwise hydrophobically lined cavity. Exit of this ion from the channel is enabled by a transient increase in the diameter of the intracellular mouth. Such "breathing" motions may form the molecular basis of channel gating.     

Molecular dynamics simulations of gramicidin A in a lipid bilayer: from structure-function relations to force fields

Molecular dynamics simulations of membrane proteins have become a popular tool for studying their dynamic features, which are not easily accessible by experiments. Whether the force fields developed for globular proteins are adequate this purpose is an important question that is often glossed over. Here we determine the permeation properties of potassium ions in the gramicidin A channel in a lipid bilayer from free energy simulations, and compare the results to experimental data. In particular, we check the dependence of the free energy barriers ions face at the channel center on the membrane size. The results indicate that there is a serious problem with the current rigid force fields independent of the membrane size, and new, possibly polarizable, force fields need to be developed to resolve this problem.     

Ion permeation through the alpha-hemolysin channel: theoretical studies based on Brownian dynamics and Poisson-Nernst-Plank electrodiffusion theory

Identification of the molecular interaction governing ion conduction through biological pores is one of the most important goals of modern electrophysiology. Grand canonical Monte Carlo Brownian dynamics (GCMC/BD) and three-dimensional Poisson-Nernst-Plank (3d-PNP) electrodiffusion algorithms offer powerful and general approaches to study of ion permeation through wide molecular pores. A detailed analysis of ion flows through the staphylococcal alpha-hemolysin channel based on series of simulations at different concentrations and transmembrane potentials is presented. The position-dependent diffusion coefficient is approximated on the basis of a hydrodynamic model. The channel conductance calculated by GCMC/BD is approximately 10% higher than (electrophysiologically measured) experimental values, whereas results from 3d-PNP are always 30-50% larger. Both methods are able to capture all important electrostatic interactions in equilibrium conditions. The asymmetric conductance upon the polarity of the transmembrane potential observed experimentally is reproduced by GCMC/BD and 3d-PNP. The separation of geometrical and energetic influence of the channel on ion conduction reveals that such asymmetries arise from the permanent charge distribution inside the pore. The major determinant of the asymmetry is unbalanced charge in the triad of polar residues D127, D128, and K131. The GCMC/BD or 3d-PNP calculations reproduce also experimental reversal potentials and permeability rations in asymmetric ionic solutions. The weak anionic selectivity of the channel results from the presence of the salt bridge between E111 and K147 in the constriction zone. The calculations also reproduce the experimentally derived dependence of the reversible potential to the direction of the salt gradient. The origin of such effect arises from the asymmetrical distribution of energetic barriers along the channel axis, which modulates the preferential ion passage in different directions.     

Valence selectivity of the gramicidin channel: a molecular dynamics free energy perturbation study.

The valence selectivity of the gramicidin channel is examined using computer simulations based on atomic models. The channel interior is modeled using a gramicidin-like periodic poly (L,D)-alanine beta-helix. Free energy perturbation calculations are performed to obtain the relative affinity of K+ and Cl- for the channel. It is observed that the interior of the gramicidin channel provides an energetically favorable interaction site for a cation but not for an anion. Relative to solvation in bulk water, the carbonyl CO oxygens can provide a favorable interaction to stabilize K+, whereas the amide NH hydrogens are much less effective in stabilizing Cl-. The results of the calculations demonstrate that, as a consequence of the structural asymmetry of the backbone charge distribution, a K+ cation can partition spontaneously from bulk water to the interior of the gramicidin channel, whereas a Cl- anion cannot.     

A combined molecular dynamics and diffusion model of single proton conduction through gramicidin

We develop a model for proton conduction through gramicidin based on the molecular dynamics simulations of Pomès and Roux (Biophys. J. 72:A246, 1997). The transport of a single proton through the gramicidin pore is described by a potential of mean force and diffusion coefficient obtained from the molecular dynamics. In addition, the model incorporates the dynamics of a defect in the hydrogen bonding structure of pore waters without an excess proton. Proton entrance and exit were not simulated by the molecular dynamics. The single proton conduction model includes a simple representation of these processes that involves three free parameters. A reasonable value can be chosen for one of these, and the other two can be optimized to yield a good fit to the proton conductance data of, Ann. N.Y. Acad. Sci. 339:8-20) for pH > or = 1.7. A sensitivity analysis shows the significance of this fit.     

Computational electrophysiology: the molecular dynamics of ion channel permeation and selectivity in atomistic detail

Presently, most simulations of ion channel function rely upon nonatomistic Brownian dynamics calculations, indirect interpretation of energy maps, or application of external electric fields. We present a computational method to directly simulate ion flux through membrane channels based on biologically realistic electrochemical gradients. In close analogy to single-channel electrophysiology, physiologically and experimentally relevant timescales are achieved. We apply our method to the bacterial channel PorB from pathogenic Neisseria meningitidis, which, during Neisserial infection, inserts into the mitochondrial membrane of target cells and elicits apoptosis by dissipating the membrane potential. We show that our method accurately predicts ion conductance and selectivity and elucidates ion conduction mechanisms in great detail. Handles for overcoming channel-related antibiotic resistance are identified.     

Studying the structure of single-component ceramide bilayers with molecular dynamics simulations using different force fields

Molecular dynamics simulations are performed on a lipid bilayer that consists of ceramide NS 24:0 in an attempt to examine several structural and physicochemical properties of the specific system. The simulations are carried out with five different force fields (OPLS, GROMOS, BERGER, CHARMM and GAFF) in order to evaluate and compare their performance in modelling lipid systems that contain ceramides. The examined properties include bilayer thickness, chain tilt, density profiles, order parameters, chain conformation, area per lipid and (intermolecular or intramolecular) hydrogen bonding between the head groups. Special focus is given to the lateral lipid arrangement. To this purpose, a method is proposed that utilises the radial distribution functions of the alkyl chains to derive quantitative information about the lateral lipid packing. In most cases, all force fields lead to similar results. For a few properties (e.g. intramolecular hydrogen bonding), there is some discrepancy between the force fields but the lack of respective experimental data does not allow an unambiguous conclusion on which force field is the most reliable.     

Ion permeation and selectivity of OmpF porin: a theoretical study based on molecular dynamics,Brownian dynamics,and continuum electrodiffusion theory

Three different theoretical approaches are used and compared to refine our understanding of ion permeation through the channel formed by OmpF porin from Escherichia coli. Those approaches are all-atom molecular dynamics (MD) in which ions, solvent, and lipids are represented explicitly, Brownian dynamics (BD) in which ions are represented explicitly, while solvent and lipids are represented as featureless dielectrics, and Poisson-Nernst-Planck (PNP) electrodiffusion theory in which both solvent and local ion concentrations are represented as a continuum. First, the ability of the different theoretical approaches in reproducing the equilibrium average ion density distribution in OmpF porin bathed by a 1M KCl symmetric salt solution is examined. Under those conditions the PNP theory is equivalent to the non-linear Poisson-Boltzmann (PB) theory. Analysis shows that all the three approaches are able to capture the important electrostatic interactions between ions and the charge distribution of the channel that govern ion permeation and selectivity in OmpF. The K(+) and Cl(-) density distributions obtained from the three approaches are very consistent with one another, which suggests that a treatment on the basis of a rigid protein and continuum dielectric solvent is valid in the case of OmpF. Interestingly, both BD and continuum electrostatics reproduce the distinct left-handed twisted ion pathways for K(+) and Cl(-) extending over the length of the pore which were observed previously in MD. Equilibrium BD simulations in the grand canonical ensemble indicate that the channel is very attractive for cations, particularly at low salt concentration. On an average there is 1.55 K(+) inside the pore in 10mM KCl. Remarkably, there is still 0.17 K(+) on average inside the pore even at a concentration as low as 1microM KCl. Secondly, non-equilibrium ion flow through OmpF is calculated using BD and PNP and compared with experimental data. The channel conductance in 0.2M and 1M KCl calculated using BD is in excellent accord with the experimental data. The calculations reproduce the experimentally well-known conductance-concentration relation and also reveal an asymmetry in the channel conductance (a larger conductance is observed under a positive transmembrane potential). Calculations of the channel conductance for three mutants (R168A, R132A, and K16A) in 1M KCl suggest that the asymmetry in the channel conductance arises mostly from the permanent charge distribution of the channel rather than the shape of the pore itself. Lastly, the calculated reversal potential in a tenfold salt gradient (0.1:1M KCl) is 27.4(+/-1.3)mV (BD) and 22.1(+/-0.6)mV (PNP), in excellent accord with the experimental value of 24.3mV. Although most of the results from PNP are qualitatively reasonable, the calculated channel conductance is about 50% higher than that calculated from BD probably because of a lack of some dynamical ion-ion correlations.     

Membrane dipole potential modulates proton conductance through gramicidin channel: movement of negative ionic defects inside the channel.

The effect of membrane dipole potential on gramicidin channel activity in bilayer lipid membranes (BLMs) was studied. Remarkably, it appeared that proton conductance of gramicidin A (gA) channels responded to modulation of the dipole potential oppositely as compared with gA alkali metal cation conductance. In particular, the addition of phloretin, known to reduce the membrane dipole potential, resulted in a decrease in gA proton conductance, on one hand, and an increase in gA alkali metal conductance, on the other hand, whereas 6-ketocholestanol, the agent raising the membrane dipole potential, provoked an increase in gA proton conductance as opposed to a decrease in the alkali metal cation conductance. The peculiarity of the 6-ketocholestanol effect consisted in its dependence on the H(+) concentration. The experiments with the impermeant dipolar compound, phloridzin, showed that the response of proton transport through gramicidin channels to varying the membrane dipole potential did not change qualitatively if the dipole potential of only one monolayer or both monolayers of the BLM was altered. In contrast to gA proton conductance, the single-channel lifetime changed similarly with varying the membrane dipole potential, regardless of the kind of permeant cations (protons or potassium ions). The results of this study could be tentatively accounted for by an assumption that one of the rate-limiting steps of proton conduction through gramicidin channels represents, in fact, movement of negatively charged species (negative ionic defects) across a membrane.     

Ion solvation and structural stability in a sodium channel investigated by molecular dynamics calculations

The stability and ion binding properties of the homo-tetrameric pore domain of a prokaryotic, voltage-gated sodium channel are studied by extensive all-atom molecular dynamics simulations, with the channel protein being embedded in a fully hydrated lipid bilayer. It is found that Na(+) ion presents in a mostly hydrated state inside the wide pore of the selectivity filter of the sodium channel, in sharp contrast to the nearly fully dehydrated state for K(+) ions in potassium channels. Our results also indicate that Na(+) ions make contact with only one or two out of the four polypeptide chains forming the selectivity filter, and surprisingly, the selectivity filter exhibits robust stability for various initial ion configurations even in the absence of ions. These findings are quite different from those in potassium channels. Furthermore, an electric field above 0.5V/nm is suggested to be able to induce Na(+) permeation through the selectivity filter.     

Effects of protein binding on a lipid bilayer containing local anesthetic articaine, and the potential of mean force calculation: a molecular dynamics simulation approach

Articaine, as a local anesthetic drug has been simulated in neutral and charged forms, and its interaction with the dimyristoylphosphatidylcholine (DMPC) lipid bilayer membrane is investigated by molecular dynamics simulation using GROMACS software. In order to obtain the optimum location of the drug molecules, as they penetrate into the membrane, umbrella sampling is applied and the free energy is calculated. The effect of protein binding to DMPC membrane on the process of drug diffusion through the membrane is considered. Five simulation systems are designed and by applying the potential of mean force, the molecular dynamics simulation on the system is performed. In light of the obtained results, the electrostatic potential, variation of lipid bilayer’s order parameter and the diffusion coefficient of drug are discussed.

Investigating various many-body force fields for their ability to predict reduction in elastic modulus due to vacancies using molecular dynamics simulations

ABSTRACT

Molecular dynamics simulations are more frequently being utilised to predict macroscale mechanical properties as a result of atomistic defects. However, the interatomic force field can significantly affect the resulting mechanical properties. While several studies exist which demonstrate the ability of various force fields to predict mechanical properties, the investigation into which is most accurate for the investigation of vacancies is limited. To obtain meaningful predictions of mechanical properties, a clear understanding of force field parameterisation is required. As such, the current study evaluates various many-body force fields to demonstrate the reduction in mechanical properties of iron and iron–chromium due to the presence of vacancies while undergoing room temperature atomistic uniaxial tension. Reduction was normalised in each case with the zero-vacancy elastic modulus, removing the need to predict an accurate nominal elastic modulus. Comparisons were made to experimental data and an empirical model from literature. It was demonstrated that accurate fitting to vacancy formation and migration energy allowed for accurate predictions. In addition, bond-order based force fields showed enhanced predictions regardless of fitting procedure. Overall, these findings highlight the need to understand capabilities and limitations of available force fields, as well as the need for enhanced parameterisation of force fields.