Comment on "Free energy simulations of single and double ion occupancy in gramicidin A" [J. Chem. Phys. 126, 105103 (2007)]

Quantitative insights into the mechanism of proton conduction and selectivity for the human voltage-gated proton channel Hv1

Human voltage-gated proton (hHv1) channels are crucial for regulating essential biological processes such as immune cell respiratory burst, sperm capacitation, and cancer cell migration. Despite the significant concentration difference between protons and other ions in physiological conditions, hHv1 demonstrates remarkable proton selectivity. Our calculations of single-proton, cation, and anion permeation free energy profiles quantitatively demonstrate that the proton selectivity of the wild-type channel originates from its strong proton affinity via the titration of the key residues D112 and D174, although the channel imposes similar kinetic blocking effects for protons compared to other ions. A two-proton knock-on model is proposed to mathematically explain the electrophysiological measurements of the pH-dependent proton conductance in the conductive state. Moreover, it is shown that the anion selectivity of the D112N mutant channel is tied to impaired proton transport and substantial anion leakage.

Accurate and Transferable Reactive Molecular Dynamics Models from Constrained Density Functional Theory

Chemical reactions constitute the central feature of many liquid, material, and biomolecular processes. Conventional molecular dynamics (MD) is inadequate for simulating chemical reactions given the fixed bonding topology of most force fields, while modeling chemical reactions using ab initio molecular dynamics is limited to shorter time and length scales given its high computational cost. As such, the multiscale reactive molecular dynamics method provides one promising alternative for simulating complex chemical systems at atomistic detail on a reactive potential energy surface. However, the parametrization of such models is a key barrier to their applicability and success. In this work, we present reactive MD models derived from constrained density functional theory that are both accurate and transferable. We illustrate the features of these models for proton dissociation reactions of amino acids in both aqueous and protein environments. Specifically, we present models for ionizable glutamate and lysine that predict accurate absolute pK_a values in water as well as their significantly shifted pK_a in staphylococcal nuclease (SNase) without any modification of the models. As one outcome of the new methodology, the simulations show that the deprotonation of ionizable residues in SNase can be closely coupled with side chain rotations, which is a concept likely generalizable to many other proteins. Furthermore, the present approach is not limited to only pK_a prediction but can enable the fully atomistic simulation of many other reactive systems along with a determination of the key aspects of the reaction mechanisms.

Polarization Effects in Water-Mediated Selective Cation Transport across a Narrow Transmembrane Channel

Despite the progress in modeling complex molecular systems of ever-increasing complexity, a quantitatively accurate computational treatment of ion permeation through narrow membrane channels remains challenging. An important factor to reach this goal is induced electronic polarization, which is likely to impact the permeation rate of small ions through narrow molecular pores. In this work, we extended the recently developed polarizable force field based on the classical Drude oscillators to assess the role of induced polarization effects on the energetics of sodium and potassium ion transport across the gramicidin A (gA) ion channel. The inclusion of induced polarization lowers barriers present in 1D potential of mean force (PMF) for cation permeation by ∼50% compared to those obtained with the additive force field. Conductance properties calculated with 1D PMFs from Drude simulations are in better agreement with experimental results. Polarization of single-file water molecules and protein atoms forming the narrow pore has a direct impact on the free-energy barriers and cation-specific solid-state NMR chemical shifts. Sensitivity analysis indicates that small changes to water-channel interactions can alter the free energy barrier for ion permeation. These results, illustrating polarization effects present in the complex electrostatic environment of the gA channel, have broad implications for revising proposed mechanisms of ion permeation and selectivity in a variety of ion channels.

Bases of Bacterial Sodium Channel Selectivity Among Organic Cations

Hille's (1971) seminal study of organic cation selectivity of eukaryotic voltage-gated sodium channels showed a sharp size cut-off for ion permeation, such that no ion possessing a methyl group was permeant. Using the prokaryotic channel, NaChBac, we found some similarity and two peculiar differences in the selectivity profiles for small polyatomic cations. First, we identified a diverse group of minimally permeant cations for wildtype NaChBac, ranging in sizes from ammonium to guanidinium and tetramethylammonium; and second, for both ammonium and hydrazinium, the charge-conserving selectivity filter mutation (E191D) yielded substantial increases in relative permeability (PX/PNa). The relative permeabilities varied inversely with relative Kd calculated from 1D Potential of Mean Force profiles (PMFs) for the single cations traversing the channel. Several of the cations bound more strongly than Na+, and hence appear to act as blockers, as well as charge carriers. Consistent with experimental observations, the E191D mutation had little impact on Na+ binding to the selectivity filter, but disrupted the binding of ammonium and hydrazinium, consequently facilitating ion permeation across the NaChBac-like filter. We concluded that for prokaryotic sodium channels, a fine balance among filter size, binding affinity, occupancy, and flexibility seems to contribute to observed functional differences.

Ligand-Binding Calculations with Metadynamics

All-atom molecular dynamics simulations can capture the dynamic degrees of freedom that characterize molecular recognition, the knowledge of which constitutes the cornerstone of rational approaches to drug design and optimization. In particular, enhanced sampling algorithms, such as metadynamics, are powerful tools to dramatically reduce the computational cost required for a mechanistic description of the binding process. Here, we describe the essential details characterizing these simulation strategies, focusing on the critical step of identifying suitable reaction coordinates, as well as on the different analysis algorithms to estimate binding affinity and residence times. We conclude with a survey of published applications that provides explicit examples of successful simulations for several targets.

The Origin of Coupled Chloride and Proton Transport in a Cl-/H+ Antiporter

The ClC family of transmembrane proteins functions throughout nature to control the transport of Cl^– ions across biological membranes. ClC-ec1 from Escherichia coli is an antiporter, coupling the transport of Cl^– and H⁺ ions in opposite directions and driven by the concentration gradients of the ions. Despite keen interest in this protein, the molecular mechanism of the Cl^–/H⁺ coupling has not been fully elucidated. Here, we have used multiscale simulation to help identify the essential mechanism of the Cl^–/H⁺ coupling. We find that the highest barrier for proton transport (PT) from the intra- to extracellular solution is attributable to a chemical reaction, the deprotonation of glutamic acid 148 (E148). This barrier is significantly reduced by the binding of Cl^– in the “central” site (Cl^–_cen), which displaces E148 and thereby facilitates its deprotonation. Conversely, in the absence of Cl^–_cen E148 favors the “down” conformation, which results in a much higher cumulative rotation and deprotonation barrier that effectively blocks PT to the extracellular solution. Thus, the rotation of E148 plays a critical role in defining the Cl^–/H⁺ coupling. As a control, we have also simulated PT in the ClC-ec1 E148A mutant to further understand the role of this residue. Replacement with a non-protonatable residue greatly increases the free energy barrier for PT from E203 to the extracellular solution, explaining the experimental result that PT in E148A is blocked whether or not Cl^–_cen is present. The results presented here suggest both how a chemical reaction can control the rate of PT and also how it can provide a mechanism for a coupling of the two ion transport processes.

Funnel metadynamics as accurate binding free-energy method

A detailed description of the events ruling ligand/protein interaction and an accurate estimation of the drug affinity to its target is of great help in speeding drug discovery strategies. We have developed a metadynamics-based approach, named funnel metadynamics, that allows the ligand to enhance the sampling of the target binding sites and its solvated states. This method leads to an efficient characterization of the binding free-energy surface and an accurate calculation of the absolute protein-ligand binding free energy. We illustrate our protocol in two systems, benzamidine/trypsin and SC-558/cyclooxygenase 2. In both cases, the X-ray conformation has been found as the lowest free-energy pose, and the computed protein-ligand binding free energy in good agreement with experiments. Furthermore, funnel metadynamics unveils important information about the binding process, such as the presence of alternative binding modes and the role of waters. The results achieved at an affordable computational cost make funnel metadynamics a valuable method for drug discovery and for dealing with a variety of problems in chemistry, physics, and material science.

Ion-controlled conformational dynamics in the outward-open transition from an occluded state of LeuT

Neurotransmitter:sodium symporter (NSS) proteins are secondary Na(+)-driven active transporters that terminate neurotransmission by substrate uptake. Despite the availability of high-resolution crystal structures of a bacterial homolog of NSSs-Leucine Transporter (LeuT)-and extensive computational and experimental structure-function studies, unanswered questions remain regarding the transport mechanisms. We used microsecond atomistic molecular-dynamics (MD) simulations and free-energy computations to reveal ion-controlled conformational dynamics of LeuT in relation to binding affinity and selectivity of the more extracellularly positioned Na(+) binding site (Na1 site). In the course of MD simulations starting from the occluded state with bound Na(+), but in the absence of substrate, we find a spontaneous transition of the extracellular vestibule of LeuT into an outward-open conformation. The outward opening is enhanced by the absence of Na1 and modulated by the protonation state of the Na1-associated Glu-290. Consistently, the Na(+) affinity for the Na1 site is inversely correlated with the extent of outward-open character and is lower than in the occluded state with bound substrate; however, the Na1 site retains its selectivity for Na(+) over K(+) in such conformational transitions. To the best of our knowledge, our findings shed new light on the Na(+)-driven transport cycle and on the symmetry in structural rearrangements for outward- and inward-open transitions.

Engineering a potent and specific blocker of voltage-gated potassium channel Kv1.3, a target for autoimmune diseases

A polypeptide toxin extracted from scorpion venom, OSK1, is modified such that its potency is drastically enhanced in blocking one class of voltage-gated potassium channels, Kv1.3, which is a pharmacological target for immunosuppressive therapy. The bound complex of Kv1.3 and OSK1 reveals that one lysine residue of the toxin is in the proximity of another lysine residue on the external vestibule of the channel, just outside of the selectivity filter. This unfavorable electrostatic interaction is eliminated by interchanging the positions of two amino acids in the toxin. The potentials of mean force of the wild-type and mutant OSK1 bound to Kv1.1-Kv1.3 channels are constructed using molecular dynamics, and the half-maximal inhibitory concentration (IC(50)) of each toxin-channel complex is computed. We show that the IC(50) values predicted for three toxins and three channels match closely with experiment. Kv1.3 is half-blocked by 0.2 pM mutant OSK1; it is >10000-fold more specific for this channel than for Kv1.1 and Kv1.2.

Gramicidin A backbone and side chain dynamics evaluated by molecular dynamics simulations and nuclear magnetic resonance experiments. I: molecular dynamics simulations

Gramicidin A (gA) channels provide an ideal system to test molecular dynamics (MD) simulations of membrane proteins. The peptide backbone lines a cation-selective pore, and due to the small channel size, the average structure and extent of fluctuations of all atoms in the peptide will influence ion permeation. This raises the question of how well molecular mechanical force fields used in MD simulations and potential of mean force (PMF) calculations can predict structure and dynamics as well as ion permeation. To address this question, we undertook a comparative study of nuclear magnetic resonance (NMR) observables predicted by fully atomistic MD simulations on a gA dimer embedded in a sodium dodecyl sulfate (SDS) micelle with measurements of the gA dimer backbone and tryptophan side chain dynamics using solution-state (15)N NMR on gA dimers in SDS micelles (Vostrikov, V. V.; Gu, H.; Ingólfsson, H. I.; Hinton, J. F.; Andersen, O. S.; Roux, B.; Koeppe, R. E., II. J. Phys. Chem. B2011, DOI 10.1021/jp200906y , accompanying article). This comparison enables us to examine the robustness of the MD simulations done using different force fields as well as their ability to predict important features of the gA channel. We find that MD is able to predict NMR observables, including the generalized order parameters (S(2)), the (15)N spin-lattice (T(1)) and spin-spin (T(2)) relaxation times, and the (1)H-(15)N nuclear Overhauser effect (NOE), with remarkable accuracy. To examine further how differences in the force fields can affect the channel conductance, we calculated the PMF for K(+) and Na(+) permeation through a gA channel in a dimyristoylphosphatidylcholine (DMPC) bilayer. In this case, we find that MD is less successful in quantitatively predicting the single-channel conductance.

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.

Making structural sense of dimerization interfaces of delta opioid receptor homodimers

Opioid receptors, like other members of the G protein-coupled receptor (GPCR) family, have been shown to associate to form dimers and/or oligomers at the plasma membrane. Whether this association is stable or transient is not known. Recent compelling evidence suggests that at least some GPCRs rapidly associate and dissociate. We have recently calculated binding affinities from free energy estimates to predict transient association between mouse delta opioid receptor (DOR) protomers at a symmetric interface involving the fourth transmembrane (TM4) helix (herein termed “4” dimer). Here we present disulfide cross-linking experiments with DOR constructs with cysteines substituted at the extracellular ends of TM4 or TM5 that confirm the formation of DOR complexes involving these helices. Our results are consistent with the involvement of TM4 and/or TM5 at the DOR homodimer interface, but possibly with differing association propensities. Coarse-grained (CG) well-tempered metadynamics simulations of two different dimeric arrangements of DOR involving TM4 alone or with TM5 (herein termed “4/5” dimer) in an explicit lipid−water environment confirmed the presence of two structurally and energetically similar configurations of the 4 dimer, as previously assessed by umbrella sampling calculations, and revealed a single energetic minimum of the 4/5 dimer. Additional CG umbrella sampling simulations of the 4/5 dimer indicated that the strength of association between DOR protomers varies depending on the protein region at the interface, with the 4 dimer being more stable than the 4/5 dimer.

Lessons from free energy simulations of delta-opioid receptor homodimers involving the fourth transmembrane helix

Several G protein-coupled receptors (GPCRs), including opioid receptors deltaOR, muOR, and kappaOR, have been reported to form stable dimers or oligomers in lipid bilayers and cell membranes. This notion has been recently challenged by imaging data supporting a transient nature of GPCR association. Here we use umbrella sampling reconstructed free energies of deltaOR homodimers involving the fourth transmembrane helix to predict their association constant. The results of these simulations, combined with estimates of diffusion-limited association rates, suggest a short lifetime for deltaOR homodimers in the membrane, in agreement with recent trends.

Using stochastic models calibrated from nanosecond nonequilibrium simulations to approximate mesoscale information

We demonstrate how the surrogate process approximation (SPA) method can be used to compute both the potential of mean force along a reaction coordinate and the associated diffusion coefficient using a relatively small number (10-20) of bidirectional nonequilibrium trajectories coming from a complex system. Our method provides confidence bands which take the variability of the initial configuration of the high-dimensional system, continuous nature of the work paths, and thermal fluctuations into account. Maximum-likelihood-type methods are used to estimate a stochastic differential equation (SDE) approximating the dynamics. For each observed time series, we estimate a new SDE resulting in a collection of SPA models. The physical significance of the collection of SPA models is discussed and methods for exploiting information in the population of estimated SPA models are demonstrated and suggested. Molecular dynamics simulations of potassium ion dynamics inside a gramicidin A channel are used to demonstrate the methodology, although SPA-type modeling has also proven useful in analyzing single-molecule experimental time series [J. Phys. Chem. B 113, 118 (2009)].

Computations of standard binding free energies with molecular dynamics simulations

An increasing number of studies have reported computations of the standard (absolute) binding free energy of small ligands to proteins using molecular dynamics (MD) simulations and explicit solvent molecules that are in good agreement with experiments. This encouraging progress suggests that physics-based approaches hold the promise of making important contributions to the process of drug discovery and optimization in the near future. Two types of approaches are principally used to compute binding free energies with MD simulations. The most widely known is the alchemical double decoupling method, in which the interaction of the ligand with its surroundings are progressively switched off. It is also possible to use a potential of mean force (PMF) method, in which the ligand is physically separated from the protein receptor. For both of these computational approaches, restraining potentials may be activated and released during the simulation for sampling efficiently the changes in translational, rotational, and conformational freedom of the ligand and protein upon binding. Because such restraining potentials add bias to the simulations, it is important that their effects be rigorously removed to yield a binding free energy that is properly unbiased with respect to the standard state. A review of recent results is presented, and differences in computational methods are discussed. Examples of computations with T4-lysozyme mutants, FKBP12, SH2 domain, and cytochrome P450 are discussed and compared. Remaining difficulties and challenges are highlighted.