Membrane assembly of simple helix homo-oligomers studied via molecular dynamics simulations

The assembly of simple transmembrane helix homo-oligomers is studied by combining a generalized Born implicit membrane model with replica exchange molecular dynamics simulations to sample the conformational space of various oligomerization states and the native oligomeric conformation. Our approach is applied to predict the structures of transmembrane helices of three proteins--glycophorin A, the M2 proton channel, and phospholamban--using only peptide sequence and the native oligomerization state information. In every case, the methodology reproduces native conformations that are in good agreement with available experimental structural data. Thus, our method should be useful in the prediction of native structures of transmembrane domains of other peptides. When we ignore the experimental constraint on the native oligomerization state and attempt de novo prediction of the structure and oligomerization state based only on sequence and simple energetic considerations, we identify the pentamer as the most stable oligomer for phospholamban. However, for the glycophorin A and the M2 proton channels, we tend to predict higher oligomers as more stable. Our studies demonstrate that reliable predictions of the structure of transmembrane helical oligomers can be achieved when the observed oligomerization state is imposed as a constraint, but that further efforts are needed for the de novo prediction of both structure and oligomeric state.     

Ion channel gating: insights via molecular simulations

Ion channels are gated, i.e. they can switch conformation between a closed and an open state. Molecular dynamics simulations may be used to study the conformational dynamics of ion channels and of simple channel models. Simulations on model nanopores reveal that a narrow (<4 A) hydrophobic region can form a functionally closed gate in the channel and can be opened by either a small (approximately 1 A) increase in pore radius or an increase in polarity. Modelling and simulation studies confirm the importance of hydrophobic gating in K channels, and support a model in which hinge-bending of the pore-lining M2 (or S6 in Kv channels) helices underlies channel gating. Simulations of a simple outer membrane protein, OmpA, indicate that a gate may also be formed by interactions of charged side chains within a pore, as is also the case in ClC channels.     

Conformational dynamics of the ligand-binding domain of inward rectifier K channels as revealed by molecular dynamics simulations: toward an understanding of Kir channel gating

Inward rectifier (Kir) potassium channels are characterized by two transmembrane helices per subunit, plus an intracellular C-terminal domain that controls channel gating in response to changes in concentration of various ligands. Based on the crystal structure of the tetrameric C-terminal domain of Kir3.1, it is possible to build a homology model of the ATP-binding C-terminal domain of Kir6.2. Molecular dynamics simulations have been used to probe the dynamics of Kir C-terminal domains and to explore the relationship between their dynamics and possible mechanisms of channel gating. Multiple simulations, each of 10 ns duration, have been performed for Kir3.1 (crystal structure) and Kir6.2 (homology model), in both their monomeric and tetrameric forms. The Kir6.2 simulations were performed with and without bound ATP. The results of the simulations reveal comparable conformational stability for the crystal structure and the homology model. There is some decrease in conformational flexibility when comparing the monomers with the tetramers, corresponding mainly to the subunit interfaces in the tetramer. The beta-phosphate of ATP interacts with the side chain of K185 in the Kir6.2 model and simulations. The flexibility of the Kir6.2 tetramer is not changed greatly by the presence of bound ATP, other than in two loop regions. Principal components analysis of the simulated dynamics suggests loss of symmetry in both the Kir3.1 and Kir6.2 tetramers, consistent with "dimer-of-dimers" motion of subunits in C-terminal domains of the corresponding Kir channels. This is suggestive of a gating model in which a transition between exact tetrameric symmetry and dimer-of-dimers symmetry is associated with a change in transmembrane helix packing coupled to gating of the channel. Dimer-of-dimers motion of the C-terminal domain tetramer is also supported by coarse-grained (anisotropic network model) calculations. It is of interest that loss of exact rotational symmetry has also been suggested to play a role in gating in the bacterial Kir homolog, KirBac1.1, and in the nicotinic acetylcholine receptor channel.     

Amino-acid solvation structure in transmembrane helices from molecular dynamics simulations

Understanding the solvation of amino acids in biomembranes is an important step to better explain membrane protein folding. Several experimental studies have shown that polar residues are both common and important in transmembrane segments, which means they have to be solvated in the hydrophobic membrane, at least until helices have aggregated to form integral proteins. In this work, we have used computer simulations to unravel these interactions on the atomic level, and classify intramembrane solvation properties of amino acids. Simulations have been performed for systematic mutations in poly-Leu helices, including not only each amino acid type, but also every z-position in a model helix. Interestingly, many polar or charged residues do not desolvate completely, but rather retain hydration by snorkeling or pulling in water/headgroups--even to the extent where many of them exist in a microscopic polar environment, with hydration levels corresponding well to experimental hydrophobicity scales. This suggests that even for polar/charged residues a large part of solvation cost is due to entropy, not enthalpy loss. Both hydration level and hydrogen bonding exhibit clear position-dependence. Basic side chains cause much less membrane distortion than acidic, since they are able to form hydrogen bonds with carbonyl groups instead of water or headgroups. This preference is supported by sequence statistics, where basic residues have increased relative occurrence at carbonyl z-coordinates. Snorkeling effects and N-/C-terminal orientation bias are directly observed, which significantly reduces the effective thickness of the hydrophobic core. Aromatic side chains intercalate efficiently with lipid chains (improving Trp/Tyr anchoring to the interface) and Ser/Thr residues are stabilized by hydroxyl groups sharing hydrogen bonds to backbone oxygens.     

Conformational analysis and molecular dynamics simulations of maltose

Constrained conformational energy minimizations have been used to calculate an adiabatic (Φ, ψ) potential energy surface for the disaccharide β-maltose. The inclusion of molecular flexibility in the conformational energy analysis of the disaccharide was found to significantly lower the barriers to conformational transitions, as has been observed previously for other systems. Several low energy wells were identified on the adiabatic surface which differ in energy by small amounts and with low absolute barriers separating them, indicating the possibility of a non-negligible equilibrium population distribution in each well. If such a distribution of conformations existed in the physical system, the conformation observed by NMR NOE measurements would thus be a “virtual” conformation. Molecular dynamics simulations of the motions of this molecule in vacuum were also conducted and indicate that the rate of relaxation of the molecule to the adiabatic surface may be slower than the typical timescale of conformational fluctuations. This effect is apparently due to an unphysical persistence of hydrogen bond patterns in vacuum which does not occur in aqueous solution. Trajectories undergoing transitions between wells were calculated and the effects of such conformational transitions upon the ensemble mean structure, such as might be observed in an NMR experiment, were demonstrated.     

OmpA: Gating and dynamics via molecular dynamics simulations

Outer membrane proteins (OMPs) of Gram-negative bacteria have a variety of functions including passive transport, active transport, catalysis, pathogenesis and signal transduction. Whilst the structures of ∼ 25 OMPs are currently known, there is relatively little known about their dynamics in different environments. The outer membrane protein, OmpA from Escherichia coli has been studied extensively in different environments both experimentally and computationally, and thus provides an ideal test case for the study of the dynamics and environmental interactions of outer membrane proteins. We review molecular dynamics simulations of OmpA and its homologues in a variety of different environments and discuss possible mechanisms of pore gating. The transmembrane domain of E. coli OmpA shows subtle differences in dynamics and interactions between a detergent micelle and a lipid bilayer environment. Simulations of the crystallographic unit cell reveal a micelle-like network of detergent molecules interacting with the protein monomers. Simulation and modelling studies emphasise the role of an electrostatic-switch mechanism in the pore-gating mechanism. Simulation studies have been extended to comparative models of OmpA homologues from Pseudomonas aeruginosa (OprF) and Pasteurella multocida (PmOmpA), the latter model including the periplasmic C-terminal domain.     

Parallel helix bundles and ion channels: molecular modeling via simulated annealing and restrained molecular dynamics.

A parallel bundle of transmembrane (TM) alpha-helices surrounding a central pore is present in several classes of ion channel, including the nicotinic acetylcholine receptor (nAChR). We have modeled bundles of hydrophobic and of amphipathic helices using simulated annealing via restrained molecular dynamics. Bundles of Ala20 helices, with N = 4, 5, or 6 helices/bundle were generated. For all three N values the helices formed left-handed coiled coils, with pitches ranging from 160 A (N = 4) to 240 A (N = 6). Pore radius profiles revealed constrictions at residues 3, 6, 10, 13, and 17. A left-handed coiled coil and a similar pattern of pore constrictions were observed for N = 5 bundles of Leu20. In contrast, N = 5 bundles of Ile20 formed right-handed coiled coils, reflecting loosened packing of helices containing beta-branched side chains. Bundles formed by each of two classes of amphipathic helices were examined: (a) M2a, M2b, and M2c derived from sequences of M2 helices of nAChR; and (b) (LSSLLSL)3, a synthetic channel-forming peptide. Both classes of amphipathic helix formed left-handed coiled coils. For (LSSLLSL)3 the pitch of the coil increased as N increased from 4 to 6. The M2c N = 5 helix bundle is discussed in the context of possible models of the pore domain of nAChR.     

Conformational dynamics of the alphaM3 transmembrane helix during acetylcholine receptor channel gating

Muscle acetylcholine receptors are synaptic ion channels that "gate" between closed- and open-channel conformations. We used Phi-value analysis to probe the transition state of the diliganded gating reaction with regard to residues in the M3, membrane-spanning helix of the muscle acetylcholine receptor alpha-subunit. Phi (a fraction between 1 and 0) parameterizes the extent to which a mutation changes the opening versus the closing rate constant and, for a linear reaction mechanism, the higher the Phi-value, the "earlier" the gating motion. In the upper half of alphaM3 the gating motions of all five tested residues were temporally correlated (Phi approximately 0.30) and serve to link structural changes occurring at the middle of the M2, pore-lining helix with those occurring at the interface of the extracellular and transmembrane domains. alphaM3 belongs to a complex and diverse set of synchronously moving parts that change structure relatively late in the channel-opening process. The propagation of the gating Brownian conformational cascade has a complex spatial distribution in the transmembrane domain.     

Defining the transmembrane helix of M2 protein from influenza A by molecular dynamics simulations in a lipid bilayer

Integral membrane proteins containing at least one transmembrane (TM) alpha-helix are believed to account for between 20% and 30% of most genomes. There are several algorithms that accurately predict the number and position of TM helices within a membrane protein sequence. However, these methods tend to disagree over the beginning and end residues of TM helices, posing problems for subsequent modeling and simulation studies. Molecular dynamics (MD) simulations in an explicit lipid and water environment are used to help define the TM helix of the M2 protein from influenza A virus. Based on a comparison of the results of five different secondary structure prediction algorithms, three different helix lengths (an 18mer, a 26mer, and a 34mer) were simulated. Each simulation system contained 127 POPC molecules plus approximately 3500-4700 waters, giving a total of approximately 18,000-21,000 atoms. Two simulations, each of 2 ns duration, were run for the 18mer and 26mer, and five separate simulations were run for the 34mer, using different starting models generated by restrained in vacuo MD simulations. The total simulation time amounted to 11 ns. Analysis of the time-dependent secondary structure of the TM segments was used to define the regions that adopted a stable alpha-helical conformation throughout the simulation. This analysis indicates a core TM region of approximately 20 residues (from residue 22 to residue 43) that remained in an alpha-helical conformation. Analysis of atomic density profiles suggested that the 18mer helix revealed a local perturbation of the lipid bilayer. Polar side chains on either side of this region form relatively long-lived H-bonds to lipid headgroups and water molecules.     

Kinked structures of isolated nicotinic receptor M2 helices: A molecular dynamics study

The pore-lining M2 helix of the nicotinic acetylcholine receptor exhibits a pronounced kink when the corresponding ion channel is in a closed conformation [N. Unwin (1993) Journal of Molecular Biology, Vol. 229, pp. 1101–1124]. We have performed molecular dynamics simulations of isolated 22-residue M2 helices in order to identify a possible molecular origin of this kink. In order to sample a wide range of conformational space, a simulated annealing protocol was used to generate five initial M2 helix structures, each of which was subsequently used as the basis of 300 ps MD simulations. Two helix sequences (M2α and M2δ) were studied in this manner, resulting in a total often 300 ps trajectories. Kinked helices present in the trajectories were identified and energy minimized to yield a total of five different stable kinked structures. For comparison, a similar molecular dynamics simulation of a Leu₂₃ helix yielded no stable kinked structures. In four of the five kinked helices, the kink was stabilized by H bonds between the helix backbone and polar side-chain atoms. Comparison with data from the literature on site-directed mutagenesis of M2 residues suggests that such polar side-chain to main-chain H bonds may also contribute to kinking of M2 helices in the intact channel protein. © 1994 John Wiley & Sons, Inc.     

Thermodynamic coupling between activation and inactivation gating in potassium channels revealed by free energy molecular dynamics simulations

The amount of ionic current flowing through K(+) channels is determined by the interplay between two separate time-dependent processes: activation and inactivation gating. Activation is concerned with the stimulus-dependent opening of the main intracellular gate, whereas inactivation is a spontaneous conformational transition of the selectivity filter toward a nonconductive state occurring on a variety of timescales. A recent analysis of multiple x-ray structures of open and partially open KcsA channels revealed the mechanism by which movements of the inner activation gate, formed by the inner helices from the four subunits of the pore domain, bias the conformational changes at the selectivity filter toward a nonconductive inactivated state. This analysis highlighted the important role of Phe103, a residue located along the inner helix, near the hinge position associated with the opening of the intracellular gate. In the present study, we use free energy perturbation molecular dynamics simulations (FEP/MD) to quantitatively elucidate the thermodynamic basis for the coupling between the intracellular gate and the selectivity filter. The results of the FEP/MD calculations are in good agreement with experiments, and further analysis of the repulsive, van der Waals dispersive, and electrostatic free energy contributions reveals that the energetic basis underlying the absence of inactivation in the F103A mutation in KcsA is the absence of the unfavorable steric interaction occurring with the large Ile100 side chain in a neighboring subunit when the intracellular gate is open and the selectivity filter is in a conductive conformation. Macroscopic current analysis shows that the I100A mutant indeed relieves inactivation in KcsA, but to a lesser extent than the F103A mutant.     

Transmembrane helix structure, dynamics, and interactions: multi-nanosecond molecular dynamics simulations.

To probe the fundamentals of membrane/protein interactions, all-atom multi-nanosecond molecular dynamics simulations were conducted on a single transmembrane poly(32)alanine helix in a fully solvated dimyristoyphosphatidylcholine (DMPC) bilayer. The central 12 residues, which interact only with the lipid hydrocarbon chains, maintained a very stable helical structure. Helical regions extended beyond these central 12 residues, but interactions with the lipid fatty-acyl ester linkages, the lipid headgroups, and water molecules made the helix less stable in this region. The C and N termini, exposed largely to water, existed as random coils. As a whole, the helix tilted substantially, from perpendicular to the bilayer plane (0 degree) to a 30 degrees tilt. The helix experienced a bend at its middle, and the two halves of the helix at times assumed substantially different tilts. Frequent hydrogen bonding, of up to 0.7 ns in duration, occurred between peptide and lipid molecules. This resulted in correlated translational diffusion between the helix and a few lipid molecules. Because of the large variation in lipid conformation, the lipid environment of the peptide was not well defined in terms of "annular" lipids and on average consisted of 18 lipid molecules. When compared with a "neat" bilayer without peptide, no significant difference was seen in the bilayer thickness, lipid conformations or diffusion, or headgroup orientation. However, the lipid hydrocarbon chain order parameters showed a significant decrease in order, especially in those methylene groups closest to the headgroup.     

Transitions from alpha to pi helix observed in molecular dynamics simulations of synthetic peptides

Molecular dynamics simulations were carried out for three 13-residue peptides of the form AcNH-A-A-E-X-A-E-A-H-A-A-E-K-A-CONH(2) with X = A, F, and W. All three peptides exhibited unexpected dynamical behavior, undergoing a transition from an alpha-helical to a pi-helical structure in the course of 5-ns trajectories in aqueous solution. Analysis of peptide length, accessible surface, interaction energies, hydrogen bonding, and dihedral angles was consistent with alpha --> pi transitions at 2800, 500, and 800 ps for X = A, F and W, respectively. The transitions occurred sequentially and cooperatively, propagating from the C- to the N-terminus for X = A and W and from the center toward both termini for X = F. The time scale of the overall transition ranged from 300 to 500 ps. For all three peptides the backbone structural transition was accompanied by a concerted rearrangement of the charged side chains, including a 3 A increase in the distance between carboxylate groups of Glu-3 and Glu-6. During the transition the peptide backbone hydrogen-bonding patterns were disrupted at the interface between the alpha-helical and nascent pi-helical regions, with peptide groups forming water-bridged hydrogen bonds. The peptide structures exhibited significant fluidity, with individual residues sampling alpha-, pi-, and 3(10)-helical conformations, as well as a "coil" state, without any intramolecular hydrogen bonds. The studied peptides have been designed to form alpha-helices when incorporated in novel hemoprotein model compounds, peptide-sandwiched mesohemes, which consist of two identical peptides covalently attached to an Fe(III) mesoporphyrin [Liu, D., Williamson, D. A., Kennedy, M. L., Williams, T. D., Morton, M. M., and Benson, D. R. (1999) J. Am. Chem. Soc. 121, 11798-11812]. The possibility of adopting pi-helical structures by the constituent peptides may influence the properties of the hemoprotein models.     

Conformational flexibility of aqueous monomeric and dimeric insulin: a molecular dynamics study

A series of molecular dynamics simulations have been used to investigate the nature of monomeric and dimeric insulin in aqueous solution. It is shown that in the absence of crystal contacts both monomeric and dimeric insulin have a high degree of intrinsic flexibility. Neither of the two monomer conformations of 2Zn crystalline insulin appears to be favored in solution nor is the asymmetry of the crystal dimer reduced in the absence of crystal contacts. A shift is observed in the relative positions of molecules 1 and 2 in the dimer compared with that found in the crystal, which may have consequences for the prediction of the effects of mutants in the monomer-monomer interface designed to alter the self-association properties of insulin.     

Efficient incorporation of protein flexibility and dynamics into molecular docking simulations

Flexibility and dynamics are protein characteristics that are essential for the process of molecular recognition. Conformational changes in the protein that are coupled to ligand binding are described by the biophysical models of induced fit and conformational selection. Different concepts that incorporate protein flexibility into protein-ligand docking within the context of these two models are reviewed. Several computational studies that discuss the validity and possible limitations of such approaches will be presented. Finally, different approaches that incorporate protein dynamics, e.g., configurational entropy, and solvation effects into docking will be highlighted.     

A new gating site in human aquaporin-4: Insights from molecular dynamics simulations

Aquaporin-4 (AQP4) is the predominant water channel in different organs and tissues. An alteration of its physiological functioning is responsible for several disorders of water regulation and, thus, is considered an attractive target with a promising therapeutic and diagnostic potential. Molecular dynamics (MD) simulations performed on the AQP4 tetramer embedded in a bilayer of lipid molecules allowed us to analyze the role of spontaneous fluctuations occurring inside the pore. Following the approach by Hashido et al. [Hashido M, Kidera A, Ikeguchi M (2007) Biophys J 93: 373–385], our analysis on 200 ns trajectory discloses three domains inside the pore as key elements for water permeation. Herein, we describe the gating mechanism associated with the well-known selectivity filter on the extracellular side of the pore and the crucial regulation ensured by the NPA motifs (asparagine, proline, alanine). Notably, on the cytoplasmic side, we find a putative gate formed by two residues, namely, a cysteine belonging to the loop D (C178) and a histidine from loop B (H95). We observed that the spontaneous reorientation of the imidazole ring of H95 acts as a molecular switch enabling H-bond interaction with C178. The occurrence of such local interaction seems to be responsible for the narrowing of the pore and thus of a remarkable decrease in water flux rate. Our results are in agreement with recent experimental observations and may represent a promising starting point to pave the way for the discovery of chemical modulators of AQP4 water permeability.     

A molecular dynamics study of gating in dioxolane-linked gramicidin A channels.

The gating transition of the RR and SS dioxolane ring-linked gramicidin A channels were studied with molecular dynamics simulations using a detailed atomic model. It was found that the probable reaction path, describing the transition of the ring from the exterior to the interior of the channel where it blocked the permeation pathway, involved several steps including the isomerization of the transpeptide plane dihedral angle of Val1. Reaction coordinates along this pathway were defined, and the transition rates between the stable conformers were calculated. It was found, in good accord with experimental observations, that the calculated blocking rate for the RR-linked channel was 280/s with a mean blocking time of 0.04 ms, whereas such blocking did not occur in the case of the SS-linked channel. An important observation is that the resulting lifetime for the blocked state of the RR-linked channel was in good accord with the experimental observations only when the calculations were performed in the presence of a potassium ion inside the channel.     

Conformational changes of pore helix coupled to gating of TRPV5 by protons

The transient receptor potential channel TRPV5 constitutes the apical entry pathway for transepithelial Ca2+ transport. We showed that TRPV5 was inhibited by both physiological intra- and extracellular acid pH. Inhibition of TRPV5 by internal protons was enhanced by extracellular acidification. Similarly, inhibition by external protons was enhanced by intracellular acidification. Mutation of either an extra- or an intracellular pH sensor blunted the cross-inhibition by internal and external protons. Both internal and external protons regulated the selectivity filter gate. Using the substituted cysteine accessibility method, we found that intracellular acidification of TRPV5 caused a conformational change of the pore helix consistent with clockwise rotation along its long axis. Thus, rotation of pore helix caused by internal protons facilitates closing of TRPV5 by external protons. This regulation by protons likely contributes to pathogenesis of disturbances of Ca2+ transport in many diseased states. Rotation of pore helix may be a common mechanism for cross-regulation of ion channels by extra- and intracellular signals.     

Conformational dynamics of chymotrypsin inhibitor 2 by coarse-grained simulations

A coarse-grained dynamic Monte Carlo (MC) simulation method is used to investigate the conformational dynamics of chymotrypsin inhibitor 2 (CI2). Each residue is represented therein by two interaction sites, one at the alpha-carbon and the other on the amino acid side-chain. The energy and geometry parameters extracted from databank structures are used. The calculated rms fluctuations of alpha-carbon atoms are in good agreement with crystallographic temperature factors. The two regions of the protein that pack against each other to form the main hydrophobic core exhibit negatively correlated fluctuations. The conformational dynamics could efficiently be probed by the time-delayed orientational and conformational correlation functions of the virtual bonds: the active site loop, excluding the active site bond, the turn region, and the N-terminal of the alpha-helix are relatively more mobile regions of the structure. A correlation is observed between the hydrogen/deuterium (H/D) exchange behavior and the long-time orientational and conformational autocorrelation function values for CI2. A cooperativity in the rotations of the bonds near in sequence is observed at all time windows, whereas the cooperative rotations of the bonds far along the sequence appear at long time windows; these correlations contribute to the stability of the secondary structures and the tertiary structure, respectively.     

Nanocontraction flows of short-chain polyethylene via molecular dynamics simulations

Nanocontraction flows of liquid short-chain polyethylene ([CH₂]₅₀) that were uniformly extruded by a constant-speed piston into a surrounding vacuum from a reservoir through an abrupt contraction nozzle were performed by employing molecular dynamics simulations. The extrudate exhibits a similar die swell phenomenon around the exit of the nozzle. In addition, numerous molecular chains are strongly adsorbed on the external surface of the nozzle. At high extrusion speeds, the velocity and temperature profiles in the nozzle show convex and concave parabolic curves, respectively, whereas the profiles are relatively flat at lower speeds. Near the internal boundary of the nozzle, the wall slip is inspected. Significantly, during the flow, the molecular chains undergo structural deformation, including compressed, stretched and shrunk motions. Comparisons with related experimental observations show that the simulated probability distributions of the bending and dihedral angles, and variations of the squared radius of gyration and orientations, are in reasonable agreement.