Ionotropic glutamate receptors are a family of tetrameric ion channels with functional states consisting of nonconducting, conducting, and desensitized states that are starting to become well characterized by electrophysiological and biophysical studies. However, the structure and relative energetics of these states beyond the general structure of the receptor are still not well understood. It is known that the interface between monomeric subunits of the tetramer plays a major role in distinguishing these functional states. We have used umbrella sampling and multimicrosecond molecular dynamics simulations of the GluA2 AMPA subtype glutamate receptor ligand-binding domain (LBD) dimers to characterize a natural propensity of the LBD dimers for various configurational states. Our results show a proposed desensitized conformation of the LBD dimer is a highly preferable conformation of the LBD dimer without the influence of other receptor domains or crystallographic conditions. This has been demonstrated by both free protein simulations of 5 μs duration, as well as by computed free energy difference between the active and desensitized states. At the same time, the simulations performed using the same protocols revealed that for the LBD mutant L483Y, known to lack desensitization, the postulated active state of the LBD dimer is indeed the preferred configurational state, which remained stable in the simulations. Our findings pave the path for developing more detailed hypotheses of the full receptor activation mechanism. Combined with the energetics of glutamate binding to the LBD and the energy required to open the transmembrane pore helices, our results strongly support a hypothesis that the low absolute free-energy state is the desensitized state of the intact AMPA receptor.

Accelerated molecular dynamics (aMD) is a promising sampling method to generate an ensemble of conformations and to explore the free energy landscape of proteins in explicit solvent. Its success resides in its ability to reduce barriers in the dihedral and the total potential energy space. However, aMD simulations of large proteins can generate large fluctuations of the dihedral and total potential energy with little conformational changes in the protein structure. To facilitate wider conformational sampling of large proteins in explicit solvent, we developed a direct intrasolute electrostatic interactions accelerated MD (DISEI-aMD) approach. This method aims to reduce energy barriers within rapidly changing electrostatic interactions between solute atoms at short-range distances. It also results in improved reconstruction quality of the original statistical ensemble of the system. Recently, we characterized a pH-dependent partial unfolding of diphtheria toxin translocation domain (T-domain) using microsecond long MD simulations. In this work, we focus on the study of conformational changes of a low-pH T-domain model in explicit solvent using DISEI-aMD. On the basis of the simulations of the low-pH T-domain model, we show that the proposed sampling method accelerates conformational rearrangement significantly faster than multiple standard aMD simulations and microsecond long conventional MD simulations.

This paper introduces explicit dependence of atomic polarizabilities on intermolecular interactions within the framework of a polarizable force field AMOEBA. Polarizable models used in biomolecular simulations often poorly describe molecular electrostatic induction in condensed phase, in part, due to neglect of a strong dependency of molecular electronic polarizability on intermolecular interactions at short distances. Our variable polarizability model parameters are derived from quantum chemical calculations of small clusters of atoms and molecules, and can be applied in simulations in condensed phase without additional scaling factors. The variable polarizability model is applied to simulate a ligand exchange reaction for a Mg2+ ion solvated in water. Explicit dependence of water polarizability on a distance between a water oxygen and Mg2+ is derived from in vacuum MP2 calculations of Mg2+–water dimer. The simulations yield a consistent description of the energetics of the Mg2+–water clusters of different size. Simulations also reproduce thermodynamics of ion solvation as well as kinetics of a water ligand exchange reaction. In contrast, simulations that used the additive force field or that used the constant polarizability models were not able to consistently and quantitatively describe the properties of the solvated Mg2+ ion.

Diphtheria toxin translocation (T) domain undergoes conformational changes in acidic solution and associates with the lipid membranes, followed by refolding and transmembrane insertion of two nonpolar helices. This process is an essential step in delivery of the toxic catalytic domain of the diphtheria toxin to the infected cell, yet its molecular determinants are poorly characterized and understood. Therefore, an atomistic model of the T-domain–membrane interaction is needed to help characterize factors responsible for such association. In this work, we present atomistic model structures of T-domain membrane-bound conformations and investigate structural factors responsible for T-domain affinity with the lipid bilayer in acidic solution using all-atom molecular dynamics (MD) simulations. The initial models of the protein conformations and protein–membrane association that serve as starting points in the present work were developed using atomistic simulations of partial unfolding of the T-domain in acidic solution (Kurnikov, I. V.; et al. J. Mol. Biol. 2013, 425, 2752–2764), and coarse-grained simulations of the T-domain association with the membranes of various compositions (Flores-Canales, J. C.; et al. J. Membr. Biol. 2015, 248, 529–543). In this work we present atomistic level modeling of two distinct configurations of the T-domain in association with the anionic lipid bilayer. In microsecond-long MD simulations both conformations retain their compact structure and gradually penetrate deeper into the bilayer interface. One membrane-bound conformation is stabilized by the protein contacts with the lipid hydrophobic core. The second modeled conformation is initially inserted less deeply and forms multiple contacts with the lipid at the interface (headgroup) region. Such contacts are formed by the charged and hydrophilic groups of partially unfolded terminal helixes and loops. Neutralization of the acidic residues at the membrane interface allows for deeper insertion of the protein and reorientation of the protein at the membrane interface, which corroborates that acidic residue protonation as well as presence of the anionic lipids may play a role in the membrane association and further membrane insertion of the T-domain as implicated in experiments. All simulations reported in this work were performed using AMBER force-field on Anton supercomputer. To perform these reported simulations, we developed and carefully tested a force-field for the anionic 1-palmitoyl-2-oleoyl-phosphatidyl-glycerol (POPG) lipid, compatible with the Amber 99SB force-field and stable in microsecond-long MD simulations in isothermal–isobaric ensemble.

Diphtheria toxin translocation (T) domain inserts in lipid bilayers upon acidification of the environment. Computational and experimental studies have suggested that low pH triggers a conformational change of the T-domain in solution preceding membrane binding. The refolded membrane-competent state was modeled to be compact and mostly retain globular structure. In the present work, we investigate how this refolded state interacts with membrane interfaces in the early steps of T-domain’s membrane association. Coarse-grained molecular dynamics simulations suggest two distinct membrane-bound conformations of the T-domain in the presence of bilayers composed of a mixture of zwitteronic and anionic phospholipids (POPC:POPG with a 1:3 molar ratio). Both membrane-bound conformations show a common near parallel orientation of hydrophobic helices TH8–TH9 relative to the membrane plane. The most frequently observed membrane-bound conformation is stabilized by electrostatic interactions between the N-terminal segment of the protein and the membrane interface. The second membrane-bound conformation is stabilized by hydrophobic interactions between protein residues and lipid acyl chains, which facilitate deeper protein insertion in the membrane interface. A theoretical estimate of a free energy of binding of a membrane-competent T-domain to the membrane is provided.

Protein molecules require both flexibility and rigidity for functioning. The fast and accurate prediction of protein rigidity/flexibility is one of the important problems in protein science. We have determined flexible regions for four homologous pairs from thermophilic and mesophilic organisms by two methods: the fast FoldUnfold which uses amino acid sequence and the time consuming MDFirst which uses three-dimensional structures. We demonstrate that both methods allow determining flexible regions in protein structure. For three of the four thermophile–mesophile pairs of proteins, FoldUnfold predicts practically the same flexible regions which have been found by the MD/First method. As expected, molecular dynamics simulations show that thermophilic proteins are more rigid in comparison to their mesophilic homologues. Analysis of rigid clusters and their decomposition provides new insights into protein stability. It has been found that the local networks of salt bridges and hydrogen bonds in thermophiles render their structure more stable with respect to fluctuations of individual contacts. Such network includes salt bridge triads Agr-Glu-Lys and Arg-Glu-Arg, or salt bridges (such as Arg-Glu) connected with hydrogen bonds. This ionic network connects alpha helices and rigidifies the structure. Mesophiles can be characterized by stand alone salt bridges and hydrogen bonds or small ionic clusters. Such difference in the network of salt bridges results in different flexibility of homologous proteins. Combining both approaches allows characterizing structural features in atomic detail that determine the rigidity/flexibility of a protein structure. This article is a part of a Special Issue entitled: The emerging dynamic view of proteins: Protein plasticity in allostery, evolution and self-assembly.Highlights► Knowledge based FoldUnfold predicts the same flexible regions as found by MDFirst. ► MDFirst uses atom contact fluctuations from MD to build graphs of protein structure. ► Rigid cluster hierarchies are different in thermophiles and mesophiles. ► In MD simulations thermophiles are less flexible then their mesophilic homologs. ► Salt bridges in protein core contribute to stability and rigidity of thermophiles.

A soft repulsion (SR) model of short-range interactions between mobile ions and protein atoms is introduced in the framework of continuum representation of the protein and solvent. The Poisson−Nernst−Plank (PNP) theory of ion transport through biological channels is modified to incorporate this soft wall protein model. Two sets of SR parameters are introduced. The first is parametrized for all essential amino acid residues using all atom molecular dynamic simulations; the second is a truncated Lennard-Jones potential. We have further designed an energy-based algorithm for the determination of the ion accessible volume, which is appropriate for a particular system discretization. The effects of these models of short-range interactions were tested by computing current−voltage characteristics of the α-hemolysin channel. The introduced SR potentials significantly improve prediction of channel selectivity. In addition, we studied the effect of the choice of some space-dependent diffusion coefficient distributions on the predicted current−voltage properties. We conclude that the diffusion coefficient distributions largely affect total currents and have little effect on rectifications, selectivity, or reversal potential. The PNP-SR algorithm is implemented in a new efficient parallel Poisson, Poisson−Boltzmann, and PNP equation solver, also incorporated in a graphical molecular modeling package HARLEM.

The ionotropic glutamate receptors are localized in the pre- and postsynaptic membrane of neurons in the brain. Activation by the principal excitatory neurotransmitter glutamate allows the ligand binding domain to change conformation, communicating opening of the channel for ion conduction. The free energy of the GluR2 S1S2 ligand binding domain (S1S2) closure transition was computed using a combination of thermodynamic integration and umbrella sampling modeling methods. A path that involves lowering the charge on E705 was chosen to clarify the role of this binding site residue. A continuum electrostatics approach in S1S2 is used to show E705, located in the ligand binding cleft, stabilizes the closed conformation of S1S2 via direct interactions with other protein residues, not through the ligand. In the closed conformation, in the absence of a ligand, S1S2 is somewhat more closed than what has been reported in X-ray structures. A semiopen conformation has been identified which is characterized by disruption of a single cross-cleft interaction and differs only slightly in energy from the fully closed S1S2. The fully open S1S2 conformation exhibits a wide energy well and shares structural similarity with the apo S1S2 crystal structure. Hybrid continuum electrostatics/MD calculations along the chosen closure transition pathway reveal solvation energies, and electrostatic interaction energies between two lobes of the protein increase the relative energetic difference between the open and closed conformational states. By analyzing the role of several cross-cleft contacts as well as other binding site residues, we demonstrate how S1S2 interactions facilitate formation of the closed conformation of the GluR2 ligand binding domain.

Conductive polymers from the polythiophene (PT) family have attracted interest in numerous domains, including potential applications in biosensing. Despite this, atomistic simulations of PTs have tended to use general organic force fields without well-tuned PT parameters, and there exists no optimized and well-validated PT force field that is compatible and consistent with existing biomolecular simulation suites. We present here the development of a new PT forcefield following the AMBER approach, using the program ANTECHAMBER and ab initio calculations at the HF/6-31G* level of theory to assign partial charges and parameterize the critical backbone torsion potential. The optimized geometries and force field potentials match well with both empirical data and previous investigators’ calculations. Initial testing of these parameters through a series of replica exchange simulations of two PT derivatives in aqueous and organic implicit solvents demonstrates that the parameters can match empirical expectations within the limits of an implicit solvent model. This new force field forms a framework for modeling of proposed PT-based devices and sensors, and is expected to accelerate device design and eventual deployment.

•Conformational switching is investigated by computational and experimental methods.•Microsecond MD simulations reveal conformational change upon histidine protonation.•TI is used to calculate free energy of histidine protonation.•Protonation of H257 ensures conformational changes.•Protonation of H223 protects from premature protonation of H257.pH-induced conformational switching is essential for functioning of diphtheria toxin, which undergoes a membrane insertion/translocation transition triggered by endosomal acidification as a key step of cellular entry. In order to establish the sequence of molecular rearrangements and side-chain protonation accompanying the formation of the membrane-competent state of the toxin's translocation (T) domain, we have developed and applied an integrated approach that combines multiple techniques of computational chemistry [e.g., long-microsecond-range, all-atom molecular dynamics (MD) simulations; continuum electrostatics calculations; and thermodynamic integration (TI)] with several experimental techniques of fluorescence spectroscopy. TI calculations indicate that protonation of H257 causes the greatest destabilization of the native structure (6.9 kcal/mol), which is consistent with our early mutagenesis results. Extensive equilibrium MD simulations with a combined length of over 8 μs demonstrate that histidine protonation, while not accompanied by the loss of structural compactness of the T-domain, nevertheless results in substantial molecular rearrangements characterized by the partial loss of secondary structure due to unfolding of helices TH1 and TH2 and the loss of close contact between the C- and N-terminal segments. The structural changes accompanying the formation of the membrane-competent state ensure an easier exposure of the internal hydrophobic hairpin formed by helices TH8 and TH9, in preparation for its subsequent transmembrane insertion.Download high-res image (325KB)Download full-size image

Tetrameric ligand binding domains of the family of ionotropic glutamate receptors assemble as dimers-of-dimers. Crystallographic studies of several glutamate receptor subtype isolated core-dimers suggest a single stable dimeric conformation. A binding domain dimer has not been captured in other conformations without the aid of biochemical methods to disrupt a critical dimer interface. Molecular dynamics simulations and continuum electrostatics calculations reveal that the active glutamate bound form of the ligand-binding domain found in typical crystal structures is the preferred energetic state of the isolated core-dimer in the presence of agonist glutamate. A desensitized conformational state is a higher energy ligand-bound state of the core-dimer. The resting apo conformational state is comparatively the least energetically favored conformation and does not contain a single state but a set of energetically equivalent conformational core-dimer states. We hypothesize the energetic balance of an open versus closed transmembrane region must be included to characterize the absolute energetic states of the full receptor, which in the presence of the ligand is believed to be a desensitized state.