Simulation of the packing of idealized transmembrane α-helix bundles

The aim of this study is to investigate if the packing motifs of native transmembrane helices can be produced by simulations with simple potentials and to develop a method for the rapid generation of initial candidate models for integral membrane proteins composed of bundles of transmembrane helices. Constituent residues are mapped along the helix axis in order to maintain the amino acid sequence-dependent properties of the helix. Helix packing is optimized according to a semi-empirical potential mainly composed of four components: a bilayer potential, a crossing angle potential, a helix dipole potential and a helix-helix distance potential. A Monte Carlo simulated annealing protocol is employed to optimize the helix bundle system. Necessary parameters are derived from theoretical studies and statistical analysis of experimentally determined protein structures. Preliminary testing of the method has been conducted with idealized seven Ala₂₀ helix bundles. The structures generated show a high degree of compactness. It was observed that both bacteriorhodopsin-like and δ-endotoxin-like structures are generated in seven-helix bundle simulations, within which the composition varies dependent upon the cooling rate. The simulation method has also been employed to explore the packing of N = 4 and N = 12 transmembrane helix bundles. The results suggest that seven and 12 transmembrane helix bundles resembling those observed experimentally (e.g., bacteriorhodopsin, rhodopsin and cytochrome c oxidase subunit I) may be generated by simulations using simple potentials. Received: 16 November 1998 / Revised version: 26 March 1999 / Accepted: 8 April 1999     

Conformational studies of human islet amyloid peptide using molecular dynamics and simulated annealing methods

Molecular dynamics simulations and simulated annealing in vacuum, model aqueous solution, and simulated membrane were used to analyze the conformational preferences of a segment spanning 20–29 residues of human islet amyloid polypeptide, [referred to as IAPP^H(20–29)]. Molecular dynamics simulations were conducted at 300 K on IAPP^H(20–29). The minimum energy conformers obtained in model aqueous solution and vacuum exhibited similar structures. Even in the absence of any constraints on peptide bonds, trans conformation was preferred consistently by all the peptide bonds. Analysis of the minimum energy conformers indicated that IAPP^H(20–29) showed a strong preference for turn structures in all the environments. These turn structures were stabilized by the formation of hydrogen bonds between the backbone amide and carbonyl groups. A good agreement was found between the results obtained from the molecular dynamics simulation and solid-state nmr experimental studies. © 1998 John Wiley & Sons, Inc. Biopoly 45: 9–20, 1998     

Improved molecular dynamics simulations for the determination of peptide structures

In this article a few methods or modifications proven to be useful in the conformational examination of peptides and related molecules by molecular dynamics are illustrated. The first is the explicit use of organic solvents in the simulations. For many cases such solvents are appropriate since the nmr measurements (or other experimental observations) were carried out in the same solvent. Here, the use of dimethylsulfoxide and chloroform in molecular dynamics is described, with some advantages of the use of these solvents highlighted. A constant allowing for the scaling of the nonbonded interactions of the force field, an idea previously employed in distance geometry and simulated annealing, has been implemented. The usefulness of this method is that when the nonbonded term is turned to zero, atoms can pass through each other, while the connectivity of the molecule is maintained. It will be shown that such simulations, if a sufficient driving force is present (i.e., nuclear Overhauser effects restraints), can produce the correct stereoconfiguration (i.e., chiral center) as well as configurational isomer (i.e., cis/trans isomers). Lastly, a penalty term for coupling constants directly related to the Karplus curve has been plemented into the potential energy force field. The advantages of this method over the commonly used dihedral angle restraining are discussed. In particular, it is shown that with more than one coupling constant about a dihedral angle a great reduction of the allowed conformational space is obtained. © 1993 John Wiley & Sons, Inc.     

Efficient modelling protocols for oligosaccharides: from vacuum to solvent

The determination of conformational preferences of oligosaccharides is best approached by describing their preferred conformations on potential energy surfaces as a function of the glycosidic linkage φ, ψ torsional angles. For proper molecular mechanics modelling the flexibility of the rotatable pendant groups must also be considered. The so called adiabatic maps partially mimic the flexibility within the 10 dimensional conformational space of the pendant groups of the given disaccharide. These molecular mechanics maps are considered to be the state-of-the art of the φ, ψ potential energy surface of disaccharides recently calculated. The RAMM (RAndom Molecular Mechanics) method was shown to be able to calculate such profiles automatically. Additionally, based on the continuum solvent approach, RAMM allows the calculation of the effects of solvent on conformational energy profiles. Molecular dynamics simulations are also useful tools to study the influence of solvent on conformational behaviour of oligosaccharides. The capability of the RAMM calculational protocol to locate low-energy conformers on the multidimensional potential energy hypersurfaces of disaccharides is illustrated and compared with molecular dynamics simulations with and without inclusion of the solvent. This revised version was published online in August 2006 with corrections to the Cover Date.     

Framework Models of Ion Permeation Through Membrane Channels and the Generalized King–Altman Method

A modern approach to studying the detailed dynamics of biomolecules is to simulate them on computers. Framework models have been developed to incorporate information from these simulations in order to calculate properties of the biomolecules on much longer time scales than can be achieved by the simulations. They also provide a simple way to think about the simulated dynamics. This article develops a method for the solution of framework models, which generalizes the King–Altman method of enzyme kinetics. The generalized method is used to construct solutions of two framework models which have been introduced previously, the single-particle and Grotthuss (proton conduction) models. The solution of the Grotthuss model is greatly simplified in comparison with direct integration. In addition, a new framework model is introduced, generalizing the shaking stack model of ion conduction through the potassium channel.     

Exhaustive conformational search and simulated annealing for models of lattice peptides

We consider simple lattice models for short peptide chains whose states can be exhaustively enumerated to find the lowest energy conformation. Using these exact results and numerical simulations, we compute the distributions for the mean time t_N, required to find the global minimum energy state by simulated annealing (SA), as a function of N, the number of units in the chain. On the basis of scaling arguments, the time t_N, to find the global minimum energy of longer chains, beyond the range covered by exhaustive enumeration, can be estimated. On the basis of the observed exponential increase in folding time of the standard SA algorithms, it is imperative that better algorithms be found for minimizing longer chains. © 1993 John Wiley & Sons, Inc.     

Molecular dynamics modeling of the sub-THz vibrational absorption of thioredoxin from <Emphasis Type="Italic">E. coli</Emphasis>

Sub-terahertz (THz) vibrational modes of the protein thioredoxin in a water environment were simulated using molecular dynamics (MD) in order to find the conditions needed for simulation convergence, improve the correlation between experimental and simulated absorption frequencies, and ultimately enhance the predictive capabilities of computational modeling. Thioredoxin from E. coli was used as a model molecule for protocol development and to optimize the simulation parameters. The empirically parameterized software packages Amber 8 and 10 were used in this work. Using atomic trajectories from the constant energy and volume MD simulations, thioredoxin’s sub-THz vibrational spectra and absorption coefficients were calculated in a quasi-harmonic approximation. An optimal production run length ∼100 ps was found, in agreement with experimental data on thioredoxin relaxation dynamics. At the same time, a new procedure was developed for averaging correlation matrices of atomic coordinates in MD simulations. In particular, the open source package ptraj was edited to improve a matrix-analyzing function. Averaging only six matrices gave much more consistent results, with absorption peak intensities exceeding those from the individual spectra and a rather good correlation between simulated vibrational frequencies and experimental data.     

Modeling the binding of three toxins to the voltage-gated potassium channel (Kv1.3)

The conduction properties of the voltage-gated potassium channel Kv1.3 and its modes of interaction with several polypeptide venoms are examined using Brownian dynamics simulations and molecular dynamics calculations. Employing an open-state homology model of Kv1.3, we first determine current-voltage and current-concentration curves and ascertain that simulated results accord with experimental measurements. We then investigate, using a molecular docking method and molecular dynamics simulations, the complexes formed between the Kv1.3 channel and several Kv-specific polypeptide toxins that are known to interfere with the conducting mechanisms of several classes of voltage-gated K⁺ channels. The depths of potential of mean force encountered by charybdotoxin, α-KTx3.7 (also known as OSK1) and ShK are, respectively, −19, −27, and −25 kT. The dissociation constants calculated from the profiles of potential of mean force correspond closely to the experimentally determined values. We pinpoint the residues in the toxins and the channel that are critical for the formation of the stable venom-channel complexes.     

Study of a structurally similar kappa opioid receptor agonist and antagonist pair by molecular dynamics simulations

Among the structurally similar guanidinonaltrindole (GNTI) compounds, 5′-GNTI is an antagonist while 6′-GNTI is an agonist of the κOR opioid receptor. To explore how a subtle alteration of the ligand structure influences the receptor activity, we investigated two concurrent processes: the final steps of ligand binding at the receptor binding site and the initial steps of receptor activation. To trace these early activation steps, the membranous part of the receptor was built on an inactive receptor template while the extracellular loops were built using the ab initio CABS method. We used the simulated annealing procedure for ligand docking and all-atom molecular dynamics simulations to determine the immediate changes in the structure of the ligand–receptor complex. The binding of an agonist, in contrast to an antagonist, induced the breakage of the “3–7 lock” between helices TM3 and TM7. We also observed an action of the extended rotamer toggle switch which suggests that those two switches are interdependent.     

On the accuracy of analytical potentials: comment on ‘Accurate ab initio calculation of the Ar–CF4 intermolecular potential energy surface’

In a recent study of the Ar–CF₄ intermolecular interaction potential [Shen C-C, Chang R-Y. Accurate ab initio calculation of the Ar–CF₄ intermolecular potential energy surface. Mol Sim. 2010;36:1111–1122], Shen and Chang (SC) illustrated how the use of bond functions can improve the accuracy and basis-set saturation of electronic structure calculations employing perturbation and coupled-cluster theory. SC then used these ab initio data to derive analytic potential energy functions for use in chemical dynamics simulations. We critically examine these analytic potentials and comment on their usage in such simulations. Our analysis highlights the need for care and global validation when deriving analytic potential energy functions.     

Constructing the suitable initial configuration of the membrane-protein system in molecular dynamics simulations

A method for constructing the suitable initial configuration of the membrane-protein system for molecular dynamics (MD) simulations is presented. This method could provide some hydrated initial configurations and help us to determine the best surface area of the system by contracting the surface area and comparing the optimized lowest energy of the system by energy minimization. The gramicidin A (GA) channel in;the fully hydrated dimyristoylphosphatidylcholine (DMPC) bilayer was used as our model. Three configurations with different surface areas were selected and applied for one 400 ps and two 300 ps MD simulations at constant pressure and temperature. All simulations were fairly stable without any constraints. Through analysis of the MD trajectories we found that the system with the best surface area was more stable than the other two systems, whose sizes were changed in the simulations. Further analysis of the bilayer normal length and the order parameters of the lipid alkyl tails indicates that the system with the best surface area shows some characteristics of the L_α phase, while both the smaller and the larger size systems have distinct deviations from the L_α phase that we expect. This illustrates that the correct surface area and the suitable initial configuration have an important influence on the phase of the membrane in the MD simulation. In addition, by comparing the root mean square differences of GA relative to the initial structure and interaction energy between different components of the system for all three systems, we find that the state of the DMPC bilayer has exerted a significant influence on the structure of GA. All these results demonstrate the validity of our method for constructing the initial configuration of the membrane-protein system for MD simulations. Received: 10 September 1998 / Revised version: 19 March 1999 / Accepted: 19 March 1999     

Stoichiometry in producer-grazer systems: Linking energy flow with element cycling

All organisms are composed of multiple chemical elements such as carbon, nitrogen and phosphorus. While energy flow and element cycling are two fundamental and unifying principles in ecosystem theory, population models usually ignore the latter. Such models implicitly assume chemical homogeneity of all trophic levels by concentrating on a single constituent, generally an equivalent of energy. In this paper, we examine ramifications of an explicit assumption that both producer and grazer are composed of two essential elements: carbon and phosphorous. Using stoichiometric principles, we construct a two-dimensional Lotka-Volterra type model that incorporates chemical heterogeneity of the first two trophic levels of a food chain. The analysis shows that indirect competition between two populations for phosphorus can shift predator—prey interactions from a (+, −) type to an unusual (−, −) class. This leads to complex dynamics with multiple positive equilibria, where bistability and deterministic extinction of the grazer are possible. We derive simple graphical tests for the local stability of all equilibria and show that system dynamics are confined to a bounded region. Numerical simulations supported by qualitative analysis reveal that Rosenzweig’s paradox of enrichment holds only in the part of the phase plane where the grazer is energy limited; a new phenomenon, the paradox of energy enrichment, arises in the other part, where the grazer is phosphorus limited. A bifurcation diagram shows that energy enrichment of producer—grazer systems differs radically from nutrient enrichment. Hence, expressing producer—grazer interactions in stoichiometrically realistic terms reveals qualitatively new dynamical behavior.     

Generalized-ensemble algorithms for molecular simulations of biopolymers.

In complex systems with many degrees of freedom such as peptides and proteins, there exists a huge number of local-minimum-energy states. Conventional simulations in the canonical ensemble are of little use, because they tend to get trapped in states of these energy local minima. A simulation in generalized ensemble performs a random walk in potential energy space and can overcome this difficulty. From only one simulation run, one can obtain canonical-ensemble averages of physical quantities as functions of temperature by the single-histogram and/or multiple-histogram reweighting techniques. In this article we review uses of the generalized-ensemble algorithms in biomolecular systems. Three well-known methods, namely, multicanonical algorithm, simulated tempering, and replica-exchange method, are described first. Both Monte Carlo and molecular dynamics versions of the algorithms are given. We then present three new generalized-ensemble algorithms that combine the merits of the above methods. The effectiveness of the methods for molecular simulations in the protein folding problem is tested with short peptide systems.     

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.     

A molecular dynamics study of a miRNA:mRNA interaction

In this paper we present a methodology to evaluate the binding free energy of a miRNA:mRNA complex through molecular dynamics (MD)–thermodynamic integration (TI) simulations. We applied our method to the Caenorhabditis elegans let-7 miRNA:lin-41 mRNA complex—a validated miRNA:mRNA interaction—in order to estimate the energetic stability of the structure. To make the miRNA:mRNA simulation possible and realistic, the methodology introduces specific solutions to overcome some of the general challenges of nucleic acid simulations and binding free energy computations that have been discussed widely in many previous research reports. The main features of the proposed methodology are: (1) positioning of the restraints imposed on the simulations in order to guarantee complex stability; (2) optimal sampling of the phase space to achieve satisfactory accuracy in the binding energy value; (3) determination of a suitable trade-off between computational costs and accuracy of binding free energy computation by the assessment of the scalability characteristics of the parallel simulations required for the TI. The experiments carried out demonstrate that MD simulations are a viable strategy for the study of miRNA binding characteristics, opening the way to the development of new computational target prediction methods based on three-dimensional structure information.     

Anomalous Surface Diffusion of Protons on Lipid Membranes

The cellular energy machinery depends on the presence and properties of protons at or in the vicinity of lipid membranes. To asses the energetics and mobility of a proton near a membrane, we simulated an excess proton near a solvated DMPC bilayer at 323 K, using a recently developed method to include the Grotthuss proton shuttling mechanism in classical molecular dynamics simulations. We obtained a proton surface affinity of −13.0 ± 0.5 kJ mol⁻¹. The proton interacted strongly with both lipid headgroup and linker carbonyl oxygens. Furthermore, the surface diffusion of the proton was anomalous, with a subdiffusive regime over the first few nanoseconds, followed by a superdiffusive regime. The time- and distance dependence of the proton surface diffusion coefficient within these regimes may also resolve discrepancies between previously reported diffusion coefficients. Our simulations show that the proton anomalous surface diffusion originates from restricted diffusion in two different surface-bound states, interrupted by the occasional bulk-mediated long-range surface diffusion. Although only a DMPC membrane was considered in this work, we speculate that the restrictive character of the on-surface diffusion is highly sensitive to the specific membrane conditions, which can alter the relative contributions of the surface and bulk pathways to the overall diffusion process. Finally, we discuss the implications of our findings for the energy machinery.     

Brownian dynamics theory for predicting internal and external blockages of tetraethylammonium in the KcsA potassium channel

The theory of Brownian dynamics is used to model permeation and the blocking of KcsA potassium channels by tetraethylammonium (TEA). A novel Brownian dynamics simulation algorithm is implemented that comprises two free energy profiles; one profile is seen by the potassium ions and the other by the TEA molecules whose shape is approximated by a sphere. Our simulations reveal that internally applied TEA blocks the passage of K⁺ ions by physically occluding the pore. A TEA molecule in the external reservoir encounters an attractive energy-well created by four tyrosine residues at position 82, in addition to all other attractive and repulsive forces impinging on it. Using Brownian dynamics, we investigate how deep the energy-well needs to be to reproduce the experimentally determined inhibitory constant k_i for the TEA blockade of KcsA or the mutant Shaker T449Y. The one-dimensional free energy profile obtained from molecular dynamics is first converted into a one-dimensional potential energy profile, and is then transformed into a three-dimensional free energy profile in Brownian dynamics by adding the short-range potential from the channel walls. When converted, the free energy profile calculated from molecular dynamics gives a well-depth of ∼10 kT. We systematically alter the depths of the profiles, and then use Brownian dynamics simulations to numerically determine the current versus TEA-concentration curves. We show that the sequence of binding and unbinding events of the TEA molecule to the binding pocket can be modeled by a first-order Markov process. The Brownian dynamics simulations also reveal that the probability of a TEA molecule binding to the binding pocket in KcsA potassium channels increases exponentially with TEA concentration and depends also on the applied potential and the K⁺ concentration in the simulation assembly.     

Insights into the functional role of protonation states in the HIV-1 protease-BEA369 complex: molecular dynamics simulations and free energy calculations

The molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method combined with molecular dynamics (MD) simulations were used to investigate the functional role of protonation in human immunodeficiency virus type 1 (HIV-1) protease complexed with the inhibitor BEA369. Our results demonstrate that protonation of two aspartic acids (Asp25/Asp25′) has a strong influence on the dynamics behavior of the complex, the binding free energy of BEA369, and inhibitor–residue interactions. Relative binding free energies calculated using the MM-PBSA method show that protonation of Asp25 results in the strongest binding of BEA369 to HIV-1 protease. Inhibitor–residue interactions computed by the theory of free energy decomposition also indicate that protonation of Asp25 has the most favorable effect on binding of BEA369. In addition, hydrogen-bond analysis based on the trajectories of the MD simulations shows that protonation of Asp25 strongly influences the water-mediated link of a conserved water molecule, Wat301. We expect that the results of this study will contribute significantly to binding calculations for BEA369, and to the design of high affinity inhibitors.