Diffraction-based density restraints for membrane and membrane-peptide molecular dynamics simulations

We have recently shown that current molecular dynamics (MD) atomic force fields are not yet able to produce lipid bilayer structures that agree with experimentally-determined structures within experimental errors. Because of the many advantages offered by experimentally validated simulations, we have developed a novel restraint method for membrane MD simulations that uses experimental diffraction data. The restraints, introduced into the MD force field, act upon specified groups of atoms to restrain their mean positions and widths to values determined experimentally. The method was first tested using a simple liquid argon system, and then applied to a neat dioleoylphosphatidylcholine (DOPC) bilayer at 66% relative humidity and to the same bilayer containing the peptide melittin. Application of experiment-based restraints to the transbilayer double-bond and water distributions of neat DOPC bilayers led to distributions that agreed with the experimental values. Based upon the experimental structure, the restraints improved the simulated structure in some regions while introducing larger differences in others, as might be expected from imperfect force fields. For the DOPC-melittin system, the experimental transbilayer distribution of melittin was used as a restraint. The addition of the peptide caused perturbations of the simulated bilayer structure, but which were larger than observed experimentally. The melittin distribution of the simulation could be fit accurately to a Gaussian with parameters close to the observed ones, indicating that the restraints can be used to produce an ensemble of membrane-bound peptide conformations that are consistent with experiments. Such ensembles pave the way for understanding peptide-bilayer interactions at the atomic level.     

Higher order multipole moments for molecular dynamics simulations

In conventional force fields, the electrostatic potential is represented by atom-centred point charges. This choice is in principle arbitrary, but technically convenient. Point charges can be understood as the first term of multipole expansions, which converge with an increasing number of terms towards the accurate representation of the molecular potential given by the electron density distribution. The use of multipole expansions can therefore improve the force field accuracy. Technically, the implementation of atomic multipoles is more involved than the use of point charges. Important points to consider are the orientation of the multipole moments during the trajectory, conformational dependence of the atomic moments and stability of the simulations which are discussed here.

Combined use of homo- and heteronuclear coupling constants as restraints in molecular dynamics simulations.

A penalty function for scalar coupling constants has been applied in molecular dynamics simulations as an experimental constraint. The function is based on the difference between the coupling constant calculated from the dihedral angle and the experimentally measured coupling constant. The method is illustrated on a model cyclic pentapeptide for which 3JHN-H alpha and 3JHN-C beta, both about the phi backbone dihedral angle, have been measured. The function is efficient in producing structures consistent with the scalar couplings, but removed from the conformation observed in solution. This arises from the lack of J restraints for the psi dihedral angle. Simulations with both nuclear Overhauser effect (NOE) and J-coupling restraints illustrates small but significant differences from simulations using only NOEs.     

On using time-averaging restraints in molecular dynamics simulation

Introducing experimental values as restraints into molecular dynamics (MD) simulations to bias the values of particular molecular properties, such as nuclear Overhauser effect intensities or distances, 3J coupling constants, chemical shifts or crystallographic structure factors, towards experimental values is a widely used structure refinement method. To account for the averaging of experimentally derived quantities inherent in the experimental techniques, time-averaging restraining methods may be used. In the case of structure refinement using 3J coupling constants from NMR experiments, time-averaging methods previously proposed can suffer from large artificially induced structural fluctuations. A modified time-averaged restraining potential energy function is proposed which overcomes this problem. The different possible approaches are compared using stochastic dynamics simulations of antamanide, a cyclic peptide of ten residues.     

Conformational regions of Boc-Ala-Aib-Ala-OMe. Sampling with molecular dynamics simulations using time averaging of distance restraints.

The method of time averaging of distance restraints in molecular dynamics simulations is applied to Boc-Ala-Aib-Ala-OMe in order to demonstrate the improved sampling properties of this method compared to conventional distance restraining. Two conformational regions, beta-turn type II and gamma-turn, are seen during MD runs at a simulation temperature of 500 K, while in simulations with conventional distance restraining, no conformational transitions could be observed for temperatures up to 1000 K.     

Electrostatic screening in molecular dynamics simulations.

The screened Coulombic potential has been shown to describe satisfactorily equilibrium properties like pK shifts, the effects of charged groups on redox potentials and binding constants of metal ions. To test how well the screening of the electrostatic potential describes the dynamical trajectory of a macromolecular system, a series of comparative simulations have been carried out on a protein system which explicitly included water molecules and a system in vacuo. For the system without solvent the results of using (i) the standard potential form were compared with results of (ii) the potential where the Coulomb term was modified by the inclusion of a distance dependent dielectric, epsilon (r), to model the screening effect of bulk water, and (iii) standard potential modified by reducing the charge on ionized residue side chains. All molecular dynamics simulations have been carried out on bovine pancreatic trypsin inhibitor. Comparisons between the resulting trajectories, averaged structures, hydrogen bonding patterns and properties such as solvent accessible surface area and radius of gyration are described. The results show that the dynamical behaviour of the protein calculated with a screened electrostatic term compares more favourably with the time-dependent structural changes of the full system with explicitly included water than the standard vacuum simulation.     

Membrane proteins: molecular dynamics simulations

Molecular dynamics simulations of membrane proteins are making rapid progress, because of new high-resolution structures, advances in computer hardware and atomistic simulation algorithms, and the recent introduction of coarse-grained models for membranes and proteins. In addition to several large ion channel simulations, recent studies have explored how individual amino acids interact with the bilayer or snorkel/anchor to the headgroup region, and it has been possible to calculate water/membrane partition free energies. This has resulted in a view of bilayers as being adaptive rather than purely hydrophobic solvents, with important implications, for example, for interaction between lipids and arginines in the charged S4 helix of voltage-gated ion channels. However, several studies indicate that the typical current simulations fall short of exhaustive sampling, and that even simple protein-membrane interactions require at least ca. 1mus to fully sample their dynamics. One new way this is being addressed is coarse-grained models that enable mesoscopic simulations on multi-mus scale. These have been used to model interactions, self-assembly and membrane perturbations induced by proteins. While they cannot replace all-atom simulations, they are a potentially useful technique for initial insertion, placement, and low-resolution refinement.     

Changes in phosphatidylcholine headgroup tilt and water order induced by monovalent salts: molecular dynamics simulations

The association between monovalent salts and neutral lipid bilayers is known to influence global bilayer structural properties such as headgroup conformational fluctuations and the dipole potential. The local influence of the ions, however, has been unknown due to limited structural resolution of experimental methods. Molecular dynamics simulations are used here to elucidate local structural rearrangements upon association of a series of monovalent Na(+) salts to a palmitoyl-oleoyl-phosphatidylcholine bilayer. We observe association of all ion types in the interfacial region. Larger anions, which are meant to rationalize data regarding a Hofmeister series of anions, bind more deeply within the bilayer than either Cl(-) or Na(+). Although the simulations are able to reproduce experimentally measured quantities, the analysis is focused on local properties currently invisible to experiments, which may be critical to biological systems. As such, for all ion types, including Cl(-), we show local ion-induced perturbations to headgroup tilt, the extent and direction of which is sensitive to ion charge and size. Additionally, we report salt-induced ordering of the water well beyond the interfacial region, which may be significant in terms of hydration repulsion between stacked bilayers.     

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.     

Study of the electrostatics treatment in molecular dynamics simulations

This article considers the treatment of long-range interactions in molecular dynamics simulations. We investigate the effects of using different cutoff distances, constant versus distance-dependent dielectric, and different smoothing methods. In contrast to findings of earlier studies, we find that increasing the cutoff over 8 A does not significantly improve the accuracy (Arnold and Ornstein, Proteins 1994;18:19-33), and using a distance-dependent dielectric instead of a constant dielectric also does not improve accuracy (Guenot and Kollman, Protein Sci 1992;1:1185-1205). This might depend on differences in simulation protocols or force fields, or both, because we use the CHARMM22 force field with stochastic boundary conditions, whereas earlier studies used other protocols and energy functions. We also note that the stability of the simulations is highly dependent on the starting structure, showing that accurate molecular simulations not only depend on a realistic simulation protocol but also on correct initial conditions.     

Nano breathers and molecular dynamics simulations in hydrogen-bonded chains

Non-linear localization phenomena in biological lattices have attracted a steadily growing interest and their existence has been predicted in a wide range of physical settings. We investigate the non-linear proton dynamics of a hydrogen-bonded chain in a semi-classical limit using the coherent state method combined with a Holstein–Primakoff bosonic representation. We demonstrate that even a weak inherent discreteness in the hydrogen-bonded (HB) chain may drastically modify the dynamics of the non-linear system, leading to instabilities that have no analog in the continuum limit. We suggest a possible localization mechanism of polarization oscillations of protons in a hydrogen-bonded chain through modulational instability analysis. This mechanism arises due to the neighboring proton–proton interaction and coherent tunneling of protons along hydrogen bonds and/or around heavy atoms. We present a detailed analysis of modulational instability, and highlight the role of the interaction strength of neighboring protons in the process of bioenergy localization. We perform molecular dynamics simulations and demonstrate the existence of nanoscale discrete breather (DB) modes in the hydrogen-bonded chain. These highly localized and long-lived non-linear breather modes may play a functional role in targeted energy transfer in biological systems.     

Active site modeling in copper azurin molecular dynamics simulations

Active site modeling in molecular dynamics simulations is investigated for the reduced state of copper azurin. Five simulation runs (5 ns each) were performed at room temperature to study the consequences of a mixed electrostatic/constrained modeling for the coordination between the metal and the polypeptide chain, using for the ligand residues a set of charges that is modified with respect to the apo form of the protein by the presence of the copper ion.The results show that the different charge values do not lead to relevant effects on the geometry of the active site of the protein, as long as bond distance constraints are used for all the five ligand atoms. The distance constraint on the O atom of Gly45 can be removed without altering the active site geometry. The coordination between Cu and the other axial ligand Met121 is outlined as being flexible. Differences are found between the bonds of the copper ion with the two apparently equivalent N¹ atoms of His46 and His117.The overall findings are discussed in connection with the issue of determining a model for the active site of azurin suitable to be used in molecular dynamics simulations under unfolding conditions. Figure Model of azurin active site. Copper ligand residues are cut off at C position except Gly45, for which the portion of backbone connecting it to His46 is shown. Only polar H atoms are shown. All atoms are in standard colors (Cu in violet), and the five ligands are labeled     

Electrophoresis in protein crystal: nonequilibrium molecular dynamics simulations

Electrophoresis of a mixture of NaCl and CaCl₂ in a lysozyme crystal is investigated using nonequilibrium molecular dynamics (MD) simulations. Upon exposure to an electric field, the stability of lysozyme is found to decrease slightly. This finding is demonstrated by increases in the root mean-square deviations of the heavy atoms of lysozyme, in the solvent-accessible surface area of hydrophobic residues, and in the number of hydrogen bonds between lysozyme and water. The solvent-accessible surface area of hydrophilic residues changes marginally, and the number of hydrogen bonds between lysozyme molecules decreases. Water molecules tend to align preferentially parallel to the electric field, and the dipole moment along the pore axis increases linearly with increasing field strength. Two pronounced layered structures are observed for Na⁺ and Ca²⁺ in the vicinity of protein surface, but only one enriched layer is observed for Cl⁻. The number distributions of all ions are nearly independent of the electric field. The water coordination numbers of all ions are smaller in the crystal than in aqueous bulk solution; however, the reverse is found for the Cl⁻ coordination numbers of cations. Both the water and the Cl⁻ coordination numbers are insensitive to the electric field. Ion diffusivities in the crystal are ∼2 orders of magnitude smaller than those in aqueous bulk solution. The drift velocities of ions increase proportionally to the electric field, particularly at high strengths, and depend on ionic charge and coordination with oppositely charged ions. Electrical current exhibits a linear relationship with the field strength. The zero-field electrical conductivity is estimated to be 0.56 S/m, which is very close to 0.61 S/m as predicted by the Nernst-Einstein equation.     

Parametrisation of time-averaged distance restraints in MD simulations

Summary Time-averaged restraints in molecular dynamics simulations offer a means to account for the averaging that is implicit in NMR spectroscopic data. We present a systematic investigation of the parameters which characterise time-averaged distance restraints. Using previously published data for a small protein, chymotrypsin inhibitor 2, we identify conditions which can lead to undesirable heating or which grossly distort the dynamics of the system.Abbreviations NOE nuclear Overhauser effect - MD molecular dynamics - CI-2 chymotrypsin inhibitor 2     

Sequence-dependent DNA deformability studied using molecular dynamics simulations

Proteins recognize specific DNA sequences not only through direct contact between amino acids and bases, but also indirectly based on the sequence-dependent conformation and deformability of the DNA (indirect readout). We used molecular dynamics simulations to analyze the sequence-dependent DNA conformations of all 136 possible tetrameric sequences sandwiched between CGCG sequences. The deformability of dimeric steps obtained by the simulations is consistent with that by the crystal structures. The simulation results further showed that the conformation and deformability of the tetramers can highly depend on the flanking base pairs. The conformations of xATx tetramers show the most rigidity and are not affected by the flanking base pairs and the xYRx show by contrast the greatest flexibility and change their conformations depending on the base pairs at both ends, suggesting tetramers with the same central dimer can show different deformabilities. These results suggest that analysis of dimeric steps alone may overlook some conformational features of DNA and provide insight into the mechanism of indirect readout during protein–DNA recognition. Moreover, the sequence dependence of DNA conformation and deformability may be used to estimate the contribution of indirect readout to the specificity of protein–DNA recognition as well as nucleosome positioning and large-scale behavior of nucleic acids.     

Application of molecular dynamics simulations to spin-labeled oligonucleotides

The EPR study of spin labeled macromolecules has provided insight into structural and dynamical properties of DNA, proteins, and related systems. While spin labeling has been useful, it is experimentally difficult to determine if the spin label significantly alters the conformation of the macromolecule to which it is attached. Molecular modeling has proven to be a powerful tool for studying structure and dynamics of biologically important molecules. Here, we have conducted molecular dynamics (MD) studies of spin labeled oligonucleotides (ONs) bearing a five (5sp) or six (6sp) membered ring nitroxide, and the corresponding unmodified ON using the suite of programs contained in Amber 5.0 with the Cornell et al. 94 force field (Cornell, W. D., Cieplak, P., Bayly, C. I., Gould, I. R., Merz, Jr., K. M. Ferguson, D. M., Spellmeyer, D. C., Fox, T., Caldwell, J. W., and Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins and Nucleic Acids. J. Am. Chem. Soc. 117, 5179-5197 (1995)). Quantum mechanical calculations employing the B3LYP method with the standard 6-31G* basis set using Gaussian98 were performed and, together with available crystallographic data for analogous nitroxides, new parameters for the nitrogen, oxygen, nitroxide alpha-carbon, and sp-hybridized carbon atoms have been developed suitable for the Cornell et al. 94 force field. MD simulations on the double-stranded (ds) spin labeled ONs, along with the corresponding unmodified analogues, have been studied over the course of 4 ns and conformational properties of all ONs are described based on the analysis of the trajectories. The spin labels were found to alter the global conformation of the ONs to which they were attached to accommodate the spin labels. The major changes include widening the major groove, decreasing helical twist, and more negative X-displacement of the base pairs. The magnitude of the effect was dependent on the specific structure of the spin label. Average and 'most representative' structures derived from the molecular dynamics simulations correlate with the experimental data on the spin labeled ONs.     

Predicting thymine dimerization yields from molecular dynamics simulations

It was recently shown that thymine dimers in the all-thymine oligonucleotide (dT)₁₈ are fully formed in <1 ps after ultraviolet excitation. The speed and low quantum yield of this reaction suggest that only a small fraction of the conformers of this structurally disordered oligonucleotide are in a position to react at the instant of photon absorption. In this work, we explore the hypothesis that conventional molecular dynamics simulations can be used to predict the yield of cyclobutane pyrimidine dimers in DNA. Conformations obtained from simulations of thymidylyl-(3′-5′)-thymidine in various cosolvents were classified as dimerizable or nondimerizable depending on the distance between the C5-C6 double bonds of the adjacent thymine bases and the torsion angle between them. The quantum yield of cyclobutane pyrimidine dimer formation was calculated as the number of dimerizable conformations divided by the total number of conformations. The experimental quantum yields measured in the different solvents were satisfactorily reproduced using physically reasonable values for the two parameters. The mean dimerizable structure computed by averaging all of the dimerizable cis-syn conformations is structurally similar to the actual cis-syn dimer. Compared to the canonical B-form TT step, the most important structural property of a dimerizable conformation is its reduced helical twist angle of 22°.     

Protein flexibility using constraints from molecular dynamics simulations

Proteins are held together in the native state by hydrophobic interactions, hydrogen bonds and interactions with the surrounding water, whose strength as well as spatial and temporal distribution affects protein flexibility and hence function. We study these effects using 10 ns molecular dynamics simulations of pure water and of two proteins, the glutamate receptor ligand binding domain and barnase. We find that most of the noncovalent interactions flicker on and off over typically nanoseconds, and so we can obtain good statistics from the molecular dynamics simulations. Based on this information, a topological network of rigid bonds corresponding to a protein structure with covalent and noncovalent bonds is constructed, with account being taken of the influence of the flickering hydrogen bonds. We define the duty cycle for the noncovalent interactions as the percentage of time a given interaction is present, which we use as an input to investigate flexibility/rigidity patterns, in the algorithm FIRST which constructs and analyses topological networks.