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

Structural, Functional and Molecular Dynamics Examination of a de novo cloned Otopetrin-like Proton Channel in crayfish

Proton channels play a crucial role in many biological functions, as they are responsible for the selective transport of protons across cell membranes. Recently, Otopetrins, a family of eukaryotic proton-selective ion channels, have attracted significant attention due to their diverse physiological roles. Despite the importance of Otopetrins, their structural and functional properties remain relatively unexplored. As a model organism, crayfish have been extensively studied to gain insights into the functioning of the nervous system. These studies cover a wide range of aspects, including the properties of individual neurons and behavioral science. However, studying the physiological systems of crayfish poses challenges for molecular research due to limited molecular sequence information available for these organisms. In the present work was identified an originally cloned mRNA, coding an Otopetrin like proton channel in the crayfish. The coded protein was modeled in silico and possible conduction mechanisms and pathways were revealed. A plasmid of the cloned mRNA was heterologously expressed in HEK293T cells. Functional experiments on transfected cells indicated that the expressed mRNA was coupled to proton conduction across the cell membrane.

Characterizing the Water Wire in the Gramicidin Channel Found by Monte Carlo Sampling Using Continuum Electrostatics and in Molecular Dynamics Trajectories with Conventional or Polarizable Force Fields

Water molecules play a key role in all biochemical processes. They help define the shape of proteins, and they are reactant or product in many reactions and are released as ligands are bound. They facilitate the transfer of protons through transmembrane proton channel, pump and transporter proteins. Continuum electrostatics (CE) force fields used by program Multiconformation CE (MCCE) capture electrostatic interactions in biomolecules with an implicit solvent, which captures the averaged solvent water equilibrium properties. Hybrid CE methods can use explicit water molecules within the protein surrounded by implicit solvent. These hybrid methods permit the study of explicit hydrogen bond networks within the protein and allow analysis of processes such as proton transfer reactions. Yet hybrid CE methods have not been rigorously tested. Here, we present an explicit treatment of water molecules in the Gramicidin A (gA) channel using MCCE and compare the resulting distributions of water molecules and key hydration features against those obtained with explicit solvent Molecular Dynamics (MD) simulations with the nonpolarizable CHARMM36 and polarizable Drude force fields. CHARMM36 leads to an aligned water wire in the channel characterized by a large absolute net water dipole moment; the MCCE and Drude analysis lead to a small net dipole moment as the water molecules change orientation within the channel. The correct orientation is not as yet known, so these calculations identify an open question.

Deciphering the Structure of the Gramicidin A Channel in the Presence of AOT Reverse Micelles in Pentane Using Molecular Dynamics Simulations

Structural studies of proteins and, in particular, integral membrane proteins (IMPs) using solution NMR spectroscopy approaches are challenging due to not only their inherent structural complexities but also the fact that they need to be solubilized in biomimetic environments (such as micelles), which enhances the slow molecular reorientation. To deal with these difficulties and increase the effective rate of molecular reorientation, the encapsulation of IMPs in the aqueous core of the reverse micelle (RM) dissolved in a low-viscosity solvent has been proven to be a viable approach. However, the effect of the reverse micelle (RM) environment on the IMP structure and function is little known. To gain insight into these aspects, this article presents a series of atomistic unconstrained molecular dynamics (MD) of a model ion channel (gramicidin A, gA) with RMs formed with anionic surfactant diacyl chain bis(2-ethylhexyl) sodium succinate (AOT) in pentane at a water-to-surfactant molar ratio (W₀) of 6. The simulations were carried out with different protocols and starting conditions for a total of 2.4 μs and were compared with other MDs used with the gA channel inserted in models of the SDS micelle or the DMPC membrane. We show here that in the presence of AOT RMs the gA dimer did not look like the "dumbbell-like" model anticipated by experiments, where the C-terminal parts of the gA are capped with two RMs and the rest of the dimer is protected from the oil solvent by the AOT acyl chains. In contrast, the MD simulations reveal that the AOT, Na⁺, and water formed two well-defined and elongated RMs attached to the C-terminal ends of the gA dimer, while the rest is in direct contact with the pentane. The initial β^6.3 secondary structure of the gA is well conserved and filled with 6-9 waters, as in SDS micelles or the DMPC membrane. Finally, the water movement inside the gA is strongly affected by the presence of RMs at each extremity, and no passage of water molecules through the gA channel is observed even after a long simulation period, whereas the opposite was found for gA in SDS and DMPC.

Umbrella sampling molecular dynamics simulations reveal concerted ion movement through G-quadruplex DNA channels

We have applied the umbrella sampling (US) method in all-atom molecular dynamics (MD) simulations to obtain potential of mean force (PMF) profiles for ion transport through three representative G-quadruplex DNA channels: [d(TG₄T)]₄, [d(G₃T₄G₄)]₂, and d[G₄(T₄G₄)₃]. The US MD results are in excellent agreement with those obtained previously with the adaptive biasing force (ABF) method. We then utilized the unique features in the US MD method to investigate multi-ion effects in [d(G₃T₄G₄)]₂ and discovered that the concerted ion movement is crucial for fully explaining the unusual experimental results on ion movement in this particular G-quadruplex system. We anticipate that these modern free-energy methods will be useful tools in evaluating ion transport properties of other G-quadruplex DNA channels.

A Polarizable Atomic Multipole-Based Force Field for Molecular Dynamics Simulations of Anionic Lipids

In all of the classical force fields, electrostatic interaction is simply treated and explicit electronic polarizability is neglected. The condensed-phase polarization, relative to the gas-phase charge distributions, is commonly accounted for in an average way by increasing the atomic charges, which remain fixed throughout simulations. Based on the lipid polarizable force field DMPC and following the same framework as Atomic Multipole Optimized Energetics for BiomoleculAr (AMOEBA) simulation, the present effort expands the force field to new anionic lipid models, in which the new lipids contain DMPG and POPS. The parameters are compatible with the AMOEBA force field, which includes water, ions, proteins, etc. The charge distribution of each atom is represented by the permanent atomic monopole, dipole and quadrupole moments, which are derived from the ab initio gas phase calculations. Many-body polarization including the inter- and intramolecular polarization is modeled in a consistent manner with distributed atomic polarizabilities. Molecular dynamics simulations of the two aqueous DMPG and POPS membrane bilayer systems, consisting of 72 lipids with water molecules, were then carried out to validate the force field parameters. Membrane width, area per lipid, volume per lipid, deuterium order parameters, electron density profile, electrostatic potential difference between the center of the bilayer and water are all calculated, and compared with limited experimental data.

Study on the Application of the Combination of TMD Simulation and Umbrella Sampling in PMF Calculation for Molecular Conformational Transitions

Free energy calculations of the potential of mean force (PMF) based on the combination of targeted molecular dynamics (TMD) simulations and umbrella samplings as a function of physical coordinates have been applied to explore the detailed pathways and the corresponding free energy profiles for the conformational transition processes of the butane molecule and the 35-residue villin headpiece subdomain (HP35). The accurate PMF profiles for describing the dihedral rotation of butane under both coordinates of dihedral rotation and root mean square deviation (RMSD) variation were obtained based on the different umbrella samplings from the same TMD simulations. The initial structures for the umbrella samplings can be conveniently selected from the TMD trajectories. For the application of this computational method in the unfolding process of the HP35 protein, the PMF calculation along with the coordinate of the radius of gyration (R_g) presents the gradual increase of free energies by about 1 kcal/mol with the energy fluctuations. The feature of conformational transition for the unfolding process of the HP35 protein shows that the spherical structure extends and the middle α-helix unfolds firstly, followed by the unfolding of other α-helices. The computational method for the PMF calculations based on the combination of TMD simulations and umbrella samplings provided a valuable strategy in investigating detailed conformational transition pathways for other allosteric processes.

Computational Investigation of the Effect of Lipid Membranes on Ion Permeation in Gramicidin A

Membrane proteins are embedded in a lipid bilayer and interact with the lipid molecules in subtle ways. This can be studied experimentally by examining the effect of different lipid bilayers on the function of membrane proteins. Understanding the causes of the functional effects of lipids is difficult to dissect experimentally but more amenable to a computational approach. Here we perform molecular dynamics simulations and free energy calculations to study the effect of two lipid types (POPC and NODS) on the conductance of the gramicidin A (gA) channel. A larger energy barrier is found for the K⁺ potential of mean force in gA embedded in POPC compared to that in NODS, which is consistent with the enhanced experimental conductance of cations in gA embedded in NODS. Further analysis of the contributions to the potential energy of K⁺ reveals that gA and water molecules in gA make similar contributions in both bilayers but there are significant differences between the two bilayers when the lipid molecules and interfacial waters are considered. It is shown that the stronger dipole moments of the POPC head groups create a thicker layer of interfacial waters with better orientation, which ultimately is responsible for the larger energy barrier in the K⁺ PMF in POPC.

Simulations of outer membrane channels and their permeability

Channels in the outer membrane of Gram-negative bacteria provide essential pathways for the controlled and unidirectional transport of ions, nutrients and metabolites into the cell. At the same time the outer membrane serves as a physical barrier for the penetration of noxious substances such as antibiotics into the bacteria. Most antibiotics have to pass through these membrane channels to either reach cytoplasmic bound targets or to further cross the hydrophobic inner membrane. Considering the pharmaceutical significance of antibiotics, understanding the functional role and mechanism of these channels is of fundamental importance in developing strategies to design new drugs with enhanced permeation abilities. Due to the biological complexity of membrane channels and experimental limitations, computer simulations have proven to be a powerful tool to investigate the structure, dynamics and interactions of membrane channels. Considerable progress has been made in computer simulations of membrane channels during the last decade. The goal of this review is to provide an overview of the computational techniques and their roles in modeling the transport across outer membrane channels. A special emphasis is put on all-atom molecular dynamics simulations employed to better understand the transport of molecules. Moreover, recent molecular simulations of ion, substrate and antibiotics translocation through membrane pores are briefly summarized. This article is part of a Special Issue entitled: Membrane Proteins edited by J.C. Gumbart and Sergei Noskov.

Theoretical Investigation of the D83V Mutation within the Myocyte-Specific Enhancer Factor-2 Beta and Its Role in Cancer

The D83V mutation in the myocyte-specific enhancer factor-2 beta (MEF2B) gene is frequently observed in lymphomas. Surprisingly, this apparent gain-of-function mutation is within a protein that is involved in the promotion of apoptosis in B cells. To investigate the oncogenic effects of this alteration and explain its predominance over other known loss-of-function mutations of MEF2B, we propose a hypothesis that this mutation influences the dynamic folding of the C-terminal loop of the N-terminal domain of MEF2B. According to our hypothesis, the mutation allows MEF2B to bind promiscuously to a wider variety of gene promoters. A large set of molecular dynamic simulations (MD) was conducted to investigate the effects of D83V mutation in silico and support the hypothesis.

Density-functional theory study of gramicidin A ion channel geometry and electronic properties

Understanding the mechanisms underlying ion channel function from the atomic-scale requires accurate ab initio modelling as well as careful experiments. Here, we present a density functional theory (DFT) study of the ion channel gramicidin A (gA), whose inner pore conducts only monovalent cations and whose conductance has been shown to depend on the side chains of the amino acids in the channel. We investigate the ground state geometry and electronic properties of the channel in vacuum, focusing on their dependence on the side chains of the amino acids. We find that the side chains affect the ground state geometry, while the electrostatic potential of the pore is independent of the side chains. This study is also in preparation for a full, linear scaling DFT study of gA in a lipid bilayer with surrounding water. We demonstrate that linear scaling DFT methods can accurately model the system with reasonable computational cost. Linear scaling DFT allows ab initio calculations with 10,000-100,000 atoms and beyond, and will be an important new tool for biomolecular simulations.

Investigation of polarization effects in the gramicidin A channel from ab initio molecular dynamics simulations

Polarization is an important component of molecular interactions and is expected to play a particularly significant role in inhomogeneous environments such as pores and interfaces. Here we investigate the effects of polarization in the gramicidin A ion channel by performing quantum mechanics/molecular mechanics molecular dynamics (MD) simulations and comparing the results with those obtained from classical MD simulations with non-polarizable force fields. We consider the dipole moments of backbone carbonyl groups and channel water molecules as well as a number of structural quantities of interest. The ab initio results show that the dipole moments of the carbonyl groups and water molecules are highly sensitive to the hydrogen bonds (H-bonds) they participate in. In the absence of a K(+) ion, water molecules in the channel are quite mobile, making the H-bond network highly dynamic. A central K(+) ion acts as an anchor for the channel waters, stabilizing the H-bond network and thereby increasing their average dipole moments. In contrast, the K(+) ion has little effect on the dipole moments of the neighboring carbonyl groups. The weakness of the ion-peptide interactions helps to explain the near diffusion-rate conductance of K(+) ions through the channel. We also address the sampling issue in relatively short ab initio MD simulations. Results obtained from a continuous 20 ps ab initio MD simulation are compared with those generated by sampling ten windows from a much longer classical MD simulation and running each window for 2 ps with ab initio MD. Both methods yield similar results for a number of quantities of interest, indicating that fluctuations are fast enough to justify the short ab initio MD simulations.

Use of Molecular Dynamics for the Refinement of an Electrostatic Model for the In Silico Design of a Polymer Antidote for the Anticoagulant Fondaparinux

Molecular dynamics (MD) simulations results are herein incorporated into an electrostatic model used to determine the structure of an effective polymer-based antidote to the anticoagulant fondaparinux. In silico data for the polymer or its cationic binding groups has not, up to now, been available, and experimental data on the structure of the polymer-fondaparinux complex is extremely limited. Consequently, the task of optimizing the polymer structure is a daunting challenge. MD simulations provided a means to gain microscopic information on the interactions of the binding groups and fondaparinux that would have otherwise been inaccessible. This was used to refine the electrostatic model and improve the quantitative model predictions of binding affinity. Once refined, the model provided guidelines to improve electrostatic forces between candidate polymers and fondaparinux in order to increase association rate constants.

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.

Engineering aspects of biological ion channels--from biosensors to computational models for permeation

This paper presents two important engineering aspects of biological ion channels-how to build sensors out of gramicidin channels and how to construct computational models for ion channel permeation. We describe our recent research in these areas, potential challenges and possible solutions.

Charge equilibration force fields for molecular dynamics simulations of lipids, bilayers, and integral membrane protein systems

With the continuing advances in computational hardware and novel force fields constructed using quantum mechanics, the outlook for non-additive force fields is promising. Our work in the past several years has demonstrated the utility of polarizable force fields, those based on the charge equilibration formalism, for a broad range of physical and biophysical systems. We have constructed and applied polarizable force fields for lipids and lipid bilayers. In this review of our recent work, we discuss the formalism we have adopted for implementing the charge equilibration (CHEQ) method for lipid molecules. We discuss the methodology, related issues, and briefly discuss results from recent applications of such force fields. Application areas include DPPC-water monolayers, potassium ion permeation free energetics in the gramicidin A bacterial channel, and free energetics of permeation of charged amino acid analogs across the water-bilayer interface. This article is part of a Special Issue entitled: Membrane protein structure and function.

Molecular dynamics simulation of hydrated DPPC monolayers using charge equilibration force fields

We present results of molecular dynamics simulations of a model DPPC-water monolayer using charge equilibration (CHEQ) force fields, which explicitly account for electronic polarization in a classical treatment of intermolecular interactions. The surface pressure, determined as the difference between the monolayer and pure water surface tensions at 323 K, is predicted to be 22.92 ±1.29 dyne/cm, just slightly below the broad range of experimental values reported for this system. The surface tension for the DPPC-water monolayer is predicted to be 42.35 ±1.16 dyne/cm, in close agreement with the experimentally determined value of 40.9 dyne/cm. This surface tension is also consistent with the value obtained from DPPC monolayer simulations using state-of-the-art nonpolarizable force fields. The current results of simulations predict a monolayer-water potential difference relative to the pure water-air interface of 0.64 ±0.02 Volts, an improved prediction compared to the fixed-charge CHARMM27 force field, yet still overestimating the experimental range of 0.3 to 0.45 Volts. As the charge equilibration model is a purely charge-based model for polarization, the current results suggest that explicitly modeled polarization effects can offer improvements in describing interfacial electrostatics in such systems.

Structural, dynamic, and electrostatic properties of fully hydrated DMPC bilayers from molecular dynamics simulations accelerated with graphical processing units (GPUs)

We present results of molecular dynamics simulations of fully hydrated DMPC bilayers performed on graphics processing units (GPUs) using current state-of-the-art non-polarizable force fields and a local GPU-enabled molecular dynamics code named FEN ZI. We treat the conditionally convergent electrostatic interaction energy exactly using the particle mesh Ewald method (PME) for solution of Poisson's Equation for the electrostatic potential under periodic boundary conditions. We discuss elements of our implementation of the PME algorithm on GPUs as well as pertinent performance issues. We proceed to show results of simulations of extended lipid bilayer systems using our program, FEN ZI. We performed simulations of DMPC bilayer systems consisting of 17,004, 68,484, and 273,936 atoms in explicit solvent. We present bilayer structural properties (atomic number densities, electron density profiles), deuterium order parameters (S(CD)), electrostatic properties (dipole potential, water dipole moments), and orientational properties of water. Predicted properties demonstrate excellent agreement with experiment and previous all-atom molecular dynamics simulations. We observe no statistically significant differences in calculated structural or electrostatic properties for different system sizes, suggesting the small bilayer simulations (less than 100 lipid molecules) provide equivalent representation of structural and electrostatic properties associated with significantly larger systems (over 1000 lipid molecules). We stress that the three system size representations will have differences in other properties such as surface capillary wave dynamics or surface tension related effects that are not probed in the current study. The latter properties are inherently dependent on system size. This contribution suggests the suitability of applying emerging GPU technologies to studies of an important class of biological environments, that of lipid bilayers and their associated integral membrane proteins. We envision that this technology will push the boundaries of fully atomic-resolution modeling of these biological systems, thus enabling unprecedented exploration of meso-scale phenomena (mechanisms, kinetics, energetics) with atomic detail at commodity hardware prices.

Comparative study of the energetics of ion permeation in Kv1.2 and KcsA potassium channels

Biological ion channels rely on a multi-ion transport mechanism for fast yet selective permeation of ions. The crystal structure of the KcsA potassium channel provided the first microscopic picture of this process. A similar mechanism is assumed to operate in all potassium channels, but the validity of this assumption has not been well investigated. Here, we examine the energetics of ion permeation in Shaker Kv1.2 and KcsA channels, which exemplify the six-transmembrane voltage-gated and two-transmembrane inward-rectifier channels. We study the feasibility of binding a third ion to the filter and the concerted motion of ions in the channel by constructing the potential of mean force for K(+) ions in various configurations. For both channels, we find that a pair of K(+) ions can move almost freely within the filter, but a relatively large free-energy barrier hinders the K(+) ion from stepping outside the filter. We discuss the effect of the CMAP dihedral energy correction that was recently incorporated into the CHARMM force field on ion permeation dynamics.

A High-Throughput Steered Molecular Dynamics Study on the Free Energy Profile of Ion Permeation through Gramicidin A

Steered molecular dynamics (SMD) simulations for the calculation of free energies are well suited for high-throughput molecular simulations on a distributed infrastructure due to the simplicity of the setup and parallel granularity of the runs. However, so far, the computational cost limited the estimation of the free energy typically over just a few pullings, thus impeding the evaluation of statistical uncertainties involved. In this work, we performed two thousand pulls for the permeation of a potassium ion in the gramicidin A pore by all-atom molecular dynamics in order to assess the bidirectional SMD protocol with a proper amount of sampling. The estimated free energy profile still shows a statistical error of several kcal/mol, while the work distributions are estimated to be non-Gaussian at pulling speeds of 10 Å/ns. We discuss the methodology and the confidence intervals in relation to increasing amounts of computed trajectories and how different permeation pathways for the potassium ion, knock-on and sideways, affect the sampling and the free energy estimation.

Energetics of double-ion occupancy in the gramicidin A channel

To understand the energetics of double-ion occupancy in gramicidin A (gA) channels, the 2D potential of mean force (PMF) is calculated for two ions at different positions along the channel axis. The cross sections of this 2D PMF are compared with available one-ion PMFs to highlight the effect of one ion on the permeation dynamics of the other. It is found that, if the first ion stays on one side in the channel, the second ion has to pass over an additional barrier to move into the outer binding site. At the same time, both outer and inner binding sites for the second ion become shallower than those in the one-ion PMF. The calculated ion-ion repulsion for a doubly occupied channel is about 2 kcal/mol, in good agreement with previous experimental estimates. The number of water molecules inside the channel and their dipole moment are calculated to interpret the energetics of double-ion occupancy. As the first ion moves into the outer binding site and then further into the channel, the oxygen atoms of the single-file water column in the channel are oriented to point toward the ion. The observed dipole moment distribution of a singly occupied channel has only one sharp peak, and the water alignment is essentially perfect once the ion is in the inner binding site. For this reason, there is an energy penalty to accommodate a second ion at the opposite end of the channel.

Simulations of a protein crystal with a high resolution X-ray structure: evaluation of force fields and water models

We use classical molecular dynamics and 16 combinations of force fields and water models to simulate a protein crystal observed by room-temperature X-ray diffraction. The high resolution of the diffraction data (0.96 Å) and the simplicity of the crystallization solution (nearly pure water) make it possible to attribute any inconsistencies between the crystal structure and our simulations to artifacts of the models rather than inadequate representation of the crystal environment or uncertainty in the experiment. All simulations were extended for 100 ns of production dynamics, permitting some long-time scale artifacts of each model to emerge. The most noticeable effect of these artifacts is a model-dependent drift in the unit cell dimensions, which can become as large as 5% in certain force fields; the underlying cause is the replacement of native crystallographic contacts with non-native ones, which can occur with heterogeneity (loss of crystallographic symmetry) in simulations with some force fields. We find that the AMBER FF99SB force field maintains a lattice structure nearest that seen in the X-ray data, and produces the most realistic atomic fluctuations (by comparison to crystallographic B-factors) of all the models tested. We find that the choice of water model has a minor effect in comparison to the choice of protein model. We also identify a number of artifacts that occur throughout all of the simulations: excessive formation of hydrogen bonds or salt bridges between polar groups and loss of hydrophobic interactions. This study is intended as a foundation for future work that will identify individual parameters in each molecular model that can be modified to improve their representations of protein structure and thermodynamics.

Synthetic chloride-selective carbon nanotubes examined by using molecular and stochastic dynamics

Synthetic channels, such as nanotubes, offer the possibility of ion-selective nanoscale pores which can broadly mimic the functions of various biological ion channels, and may one day be used as antimicrobial agents, or for treatment of cystic fibrosis. We have designed a carbon nanotube that is selectively permeable to anions. The virtual nanotubes are constructed from a hexagonal array of carbon atoms (graphene) rolled up to form a tubular structure, with an effective radius of 4.53 Å and length of 34 Å. The pore ends are terminated with polar carbonyl groups. The nanotube thus formed is embedded in a lipid bilayer and a reservoir containing ionic solutions is added at each end of the pore. The conductance properties of these synthetic channels are then examined with molecular and stochastic dynamics simulations. Profiles of the potential of mean force at 0 mM reveal that a cation moving across the pore encounters an insurmountable free energy barrier of ∼25 kT in height. In contrast, for anions, there are two energy wells of ∼12 kT near each end of the tube, separated by a central free energy barrier of 4 kT. The conductance of the pore, with symmetrical 500 mM solutions in the reservoirs, is 72 pS at 100 mV. The current saturates with an increasing ionic concentration, obeying a Michaelis-Menten relationship. The pore is normally occupied by two ions, and the rate-limiting step in conduction is the time taken for the resident ion near the exit gate to move out of the energy well.

The gramicidin channel ion permeation free-energy profile: direct and indirect effects of CHARMM force field improvements

A revised CHARMM force field for tryptophan residues is studied as well as a new grid-based correction algorithm, called CMAP, using molecular dynamics simulations of gramicidin A (1JNO) embedded in a lipid bilayer (DMPC) with 1 mol/kg NaCl or KCl saline solution. The conformational stability of the interfacial side chains is studied, which shows good stability on the 10 ns time scale. The revised force field for the tryptophan side chain produces, in the decomposition, a Na(+) PMF(Trp) profile that is consonant with the prediction from the experimental results, analyzed with rate theory by Durrant et al. (2006), but in stark contrast to the prediction of the original CHARMM force field, version 22. However, the effect is diluted in the PMF profile due to indirect effects mediated by other components of the system (polypeptide, lipid molecules, ions, and water molecules). CMAP corrections to the L-amino acids help reduce the excessive translocation barrier. Decomposition demonstrates that this effect is due to effects on the K(+) PMF(H(2)O) profile rather than on the K(+) PMF(gA) profile. The results have been confirmed to be robust using an alternative umbrella-potential method. Further force field balancing efforts (direct and indirect) are required for future studies to evaluate whether these effects give rise to predictions that are consistent with those observables extracted from real experiments.

P-SPLINES USING DERIVATIVE INFORMATION

Time series associated with single-molecule experiments and/or simulations contain a wealth of multiscale information about complex biomolecular systems. We demonstrate how a collection of Penalized-splines (P-splines) can be useful in quantitatively summarizing such data. In this work, functions estimated using P-splines are associated with stochastic differential equations (SDEs). It is shown how quantities estimated in a single SDE summarize fast-scale phenomena, whereas variation between curves associated with different SDEs partially reflects noise induced by motion evolving on a slower time scale. P-splines assist in "semiparametrically" estimating nonlinear SDEs in situations where a time-dependent external force is applied to a single-molecule system. The P-splines introduced simultaneously use function and derivative scatterplot information to refine curve estimates. We refer to the approach as the PuDI (P-splines using Derivative Information) method. It is shown how generalized least squares ideas fit seamlessly into the PuDI method. Applications demonstrating how utilizing uncertainty information/approximations along with generalized least squares techniques improve PuDI fits are presented. Although the primary application here is in estimating nonlinear SDEs, the PuDI method is applicable to situations where both unbiased function and derivative estimates are available.

Free-energy profiles for ions in the influenza M2-TMD channel

M(2) transmembrane domain channel (M(2)-TMD) permeation properties are studied using molecular dynamics simulations of M(2)-TMD (1NYJ) embedded in a lipid bilayer (DMPC) with 1 mol/kg NaCl or KCl saline solution. This study allows examination of spontaneous cation and anion entry into the selectivity filter. Three titration states of the M(2)-TMD tetramer are modeled for which the four His(37) residues, forming the selectivity filter, are net uncharged, +2 charged, or +3 charged. M(2)-TMD structural properties from our simulations are compared with the properties of other models extracted from NMR and X-ray studies. During 10 ns simulations, chloride ions occasionally occupy the positively-charged selectivity filter region, and from umbrella sampling simulations, Cl(-) has a lower free-energy barrier in the selectivity-filter region than either Na(+) or NH(4) (+), and NH(4) (+) has a lower free-energy barrier than Na(+). For Na(+) and Cl(-), the free-energy barriers are less than 5 kcal/mol, suggesting that the 1NYJ conformation would probably not be exquisitely proton selective. We also point out a rotameric configuration of Trp(41) that could fully occlude the channel.

Charge equilibration force fields for lipid environments: applications to fully hydrated DPPC bilayers and DMPC-embedded gramicidin A

Polarizable force fields for lipid and solvent environments are used for molecular dynamics simulations of a fully hydrated dipalmitoylphosphatidylcholine (DPPC) bilayer and gramicidin A (gA) dimer embedded in a dimyristoylphosphatidylcholine (DMPC) bilayer. The lipid bilayer is modeled using the CHARMM charge equilibration (CHEQ) polarizable force field for lipids and the TIP4P-FQ force field to represent solvent. For the DPPC bilayer system, results are compared to the same system simulated using the nonpolarizable CHARMM27r (C27r) force field and TIP3P water. Calculated atomic and electron density profiles, head group orientations as measured by the phosphorus-nitrogen vector orientation, and deuterium order parameters are found to be consistent with previous simulations and with experiment. The CHEQ model exhibits greater water penetration into the bilayer interior, as demonstrated by the potential of mean force calculated from the water density profile. This is a result of the variation of the water molecular dipole from 2.55 D in the bulk to 1.88 D in the interior. We discuss this finding in the context of previous studies (both simulation and experiment) that have investigated the extent of penetration of water into DPPC bilayers. We also discuss the effects of including explicit polarization on the water dipole moment variation as a function of distance from the bilayer. We show distributions of atomic charges over the course of the simulation since the CHEQ model allows the charges to fluctuate. We have calculated the interfacial dipole potential, which the CHEQ model predicts to be 0.95 V compared to 0.86 V as predicted by the C27r model. We also discuss dielectric permittivity profiles and the differences arising between the two models. We obtain bulk values of 72.77 for the CHEQ model (TIP4P-FQ water) and 91.22 for C27r (TIP3P), and values approaching unity in the membrane interior. Finally, we present results of simulations of gA embedded in a DMPC bilayer using the CHEQ model and discuss structural properties.

Mechanism and energetics of charybdotoxin unbinding from a potassium channel from molecular dynamics simulations

Ion channel-toxin complexes are ideal systems for computational studies of protein-ligand interactions, because, in most cases, the channel axis provides a natural reaction coordinate for unbinding of a ligand and a wealth of physiological data is available to check the computational results. We use a recently determined structure of a potassium channel-charybdotoxin complex in molecular dynamics simulations to investigate the mechanism and energetics of unbinding. Pairs of residues on the channel protein and charybdotoxin that are involved in the binding are identified, and their behavior is traced during umbrella-sampling simulations as charybdotoxin is moved away from the binding site. The potential of mean force for the unbinding of charybdotoxin is constructed from the umbrella sampling simulations using the weighted histogram analysis method, and barriers observed are correlated with specific breaking of interactions and influx of water molecules into the binding site. Charybdotoxin is found to undergo conformational changes as a result of the reaction coordinate choice--a nontrivial decision for larger ligands--which we explore in detail, and for which we propose solutions. Agreement between the calculated and the experimental binding energies is obtained once the energetic consequences of these conformational changes are included in the calculations.

Computational studies of gramicidin permeation: an entry way sulfonate enhances cation occupancy at entry sites

The impact on the cation-transport free-energy profile of replacing the C-terminal ethanolamine in the gramicidin A channel with a taurine residue is studied using molecular dynamics simulations of gramicidin A (1JNO) embedded in a lipid bilayer (DMPC) with 1 mol/kg NaCl saline solution. The potential of mean force for ion transport is obtained by umbrella sampling. The presence of a negatively charged sulfonate group at the entrance of the gramicidin channel affects the depth and the location of the binding sites, producing a strong attraction for the cations in the bulk. The potential of mean force by the sulfonate acting directly through electrostatics and van der Waals interactions on the test ion is highly modulated by indirect effects (i.e., sulfonate effects on other components of the system that, in turn, affect the ion free-energy profile in the channel). Because the "entry" sites are located symmetrically at both entry and exit of the channel, the deeper free-energy wells should inhibit exit. Given that the channel has increased conductance experimentally, the simulation results suggest that the channel conductance is normally entry limited.

Molecular dynamics simulations of a DMPC bilayer using nonadditive interaction models

We present a polarizable force field based on the charge-equilibration formalism for molecular dynamics simulations of phospholipid bilayers. We discuss refinement of headgroup dihedral potential parameters to reproduce ab initio conformational energies of dimethylphosphate calculated at the MP2/cc-pVTZ level of theory. We also address the refinement of electrostatic and Lennard-Jones (van der Waals) parameters to reproduce ab initio polarizabilities and water interaction energies of dimethylphosphate and tetramethylammonium. We present results of molecular dynamics simulations of a solvated dimyristoylphosphatidylcholine bilayer using this polarizable force field as well as the nonpolarizable, fixed-charge CHARMM27 and CHARMM27r force fields for comparison. Calculated atomic and electron-density profiles, deuterium order parameters, and headgroup orientations are found to be consistent with previous simulations and with experiment. Polarizable interaction models for solvent and lipid exhibit greater water penetration into the lipid interior; this is due to the variation of water molecular dipole moment from a bulk value of 2.6 Debye to a value of 1.9 Debye in the membrane interior. The reduction in the electrostatic component of the desolvation free-energy penalty allows for greater water density. The surface dipole potential predicted by the polarizable model is 0.95 V compared to the value of 0.8 V based on nonpolarizable force-field calculations. Effects of inclusion of explicit polarization are discussed in relation to water dipole moment and varying charge distributions. Dielectric permittivity profiles for polarizable and nonpolarizable interactions exhibit subtle differences arising from the nature of the individual component parameterizations; for the polarizable force field, we obtain a bulk dielectric permittivity of 79, whereas the nonpolarizable force field plateaus at 97 (the value for pure TIP3P water). In the membrane interior, both models predict unit permittivities, with the polarizable models contributing from one to two more units due to the optical dielectric (high-frequency dipole fluctuations). This contribution is a step toward the continuing development of a CHARMM (Chemistry at Harvard Molecular Mechanics) polarizable force field for simulations of biomacromolecular systems.

Potential of mean force for Syrian hamster prion epitope protein--monoclonal fab 3f4 antibody interaction studies

Simulating antigen-antibody interactions are crucial for understanding antigen-antibody associations in immunology. To shed further light on this question, we study a dissociation of the Syrian hamster prion epitope protein-fab 3f4 antibody complex structure. The stretching, that is, the distance between the center of mass of the prion epitope protein and the fab 3f4 antibody, has been studied using potential of mean force (PMF) calculations based on molecular dynamics (MD) and the implicit water model. For the complex structure, there are four important intermediates, U-shaped groove on the antibodies, and two inter-protein molecular hydrogen bonds in the stretching process. Use of our simulations may help in understanding the binding mechanics of the complex structure, and thus of significance in the design of antibodies against prion disease.

Potential of mean force calculations of ligand binding to ion channels from Jarzynski's equality and umbrella sampling

Potential of mean force (PMF) calculations provide a reliable method for determination of the absolute binding free energies for protein-ligand systems. The common method used for this purpose -- umbrella sampling with weighted histogram analysis -- is computationally very laborious, which limits its applications. Recently, a much simpler alternative for PMF calculations has become available, namely, using Jarzynski's equality in steered molecular dynamics simulations. So far, there have been a few comparisons of the two methods and mostly in simple systems that do not reflect the complexities of protein-ligand systems. Here, we use both methods to calculate the PMF for ion permeation and ligand binding to ion channels. Comparison of results indicate that Jarzynski's method suffers from relaxation problems in complex systems and would require much longer simulation times to yield reliable PMFs for protein-ligand systems.