Assessing and refining molecular dynamics simulations of proteins with nuclear magnetic resonance data

The sophistication of the force fields, algorithms and hardware used for molecular dynamics (MD) simulations of proteins is continuously increasing. No matter how advanced the methodology, however, it is essential to evaluate the appropriateness of the structures sampled in a simulation by comparison with quantitative experimental data. Solution nuclear magnetic resonance (NMR) data are particularly useful for checking the quality of protein simulations, as they provide both structural and dynamic information on a variety of temporal and spatial scales. Here, various features and implications of using NMR data to validate and bias MD simulations are outlined, including an overview of the different types of NMR data that report directly on structural properties and of relevant simulation techniques. The focus throughout is on how to properly account for conformational averaging, particularly within the context of the assumptions inherent in the relationships that link NMR data to structural properties.     

Determination of protein structures from nuclear magnetic resonance data using a restrained molecular dynamics approach: the lac repressor DNA binding domain

A procedure is described to determine from NMR data the three-dimensional structure of biomolecules. This procedure combines model building with a restrained Molecular Dynamics algorithm, in which distance information from NOEs is incorporated in the form of pseudo potentials. The method has been applied to the N-terminal DNA-binding domain or "headpiece" (amino acids 1-51) of the lac repressor from E. coli, for which no crystal structure is available. The spatial structure of the headpiece is discussed in terms of known physical and biochemical data and of its DNA binding properties.     

Similarity of force-induced unfolding of apomyoglobin to its chemical-induced unfolding: an atomistic molecular dynamics simulation approach

We have compared force-induced unfolding with traditional unfolding methods using apomyoglobin as a model protein. Using molecular dynamics simulation, we have investigated the structural stability as a function of the degree of mechanical perturbation. Both anisotropic perturbation by stretching two terminal atoms and isotropic perturbation by increasing the radius of gyration of the protein show the same key event of force-induced unfolding. Our primary results show that the native structure of apomyoglobin becomes destabilized against the mechanical perturbation as soon as the interhelical packing between the G and H helices is broken, suggesting that our simulation results share a common feature with the experimental observation that the interhelical contact is more important for the folding of apomyoglobin than the stability of individual helices. This finding is further confirmed by simulating both helix destabilizing and interhelical packing destabilizing mutants.     

Solution conformations of a trimannoside from nuclear magnetic resonance and molecular dynamics simulations

N-linked oligosaccharides often act as ligands for receptor proteins in a variety of cell recognition processes. Knowledge of the solution conformations, as well as protein-bound conformations, of these oligosaccharides is required to understand these important interactions. In this paper we present a model for the solution conformations sampled by a simple trimannoside, methyl 3, 6-di-O-(alpha-D-mannopyranosyl)-alpha-D-mannopyranoside, which contains two of the most commonly found glycosidic linkages in N-linked oligosaccharides. This model was derived from simulated annealing protocols incorporating distance restraints extracted from NOESY spectra along with torsional restraints computed from three-bond (1)H-(13)C coupling constants measured across the glycosidic bonds. The model was refined in light of unrestrained molecular dynamics simulations conducted in the presence of solvent water. The resulting model depicts a molecule undergoing conformational averaging in solution, adopting four major and two minor conformations. The four major conformations arise from a pair of two-state transitions, one each at the alpha(1-->3) and alpha(1-->6) linkages, whereas the minor conformations result from an additional transition of the alpha(1-->6) linkage. Our data also suggest that the alpha(1-->3) transition is fast and changes the molecular shape slightly, whereas the alpha(1-->6) is much slower and alters the molecular shape dramatically.     

Mapping the intrinsically disordered properties of the flexible loop domain of Bcl-2: a molecular dynamics simulation study

Most of the B-cell lymphoma-2 (Bcl-2) protein structure has been elucidated; however, the conformation of its flexible loop domain (FLD) has not yet been experimentally predicted. Its high flexibility under physiological conditions is the reason. FLD behaves as an intrinsically disordered region (IDR) and can adopt regular structures in particular conditions associated with the control of Bcl-2’s anti-apoptotic functions. In a previous contribution, we analyzed an engineered Bcl-2 construct (Bcl-2-Δ22Σ3) submitted to 25-ns MD and reported a disordered-to-helix transitions in a region of FLD (rFLD, residues 60–77). However, the conformational preferences in solution of rFLD in the nanosecond to microsecond scale were not analyzed. Herein, an average model was obtained for the native Bcl-2 protein by homology modeling and MD simulation techniques. From this, only the atomic coordinates corresponding to the rFLD were simulated for 1 μs by MD at 310 K. In concordance with previous studies, a disordered-to-helix transitions were exhibited, implying that this “interconversion of folding” in the rFLD suggest a possible set of conformations encoded in its sequence. Principal component analysis (PCA) showed that most of the conformational fluctuation of Bcl-2 is provided by rFLD. Dihedral PCA (dPCA) offered information about all the conformations of rFLD in the μs of the simulation, characterizing a dPCA-based free energy landscape of rFLD, and a conformational ensemble of fast interconverting conformations as other IDRs. Furthermore, despite the conformational heterogeneity of rFLD, the analysis of the dihedral angles (Φ, Ψ) showed that this region does not randomly explore the conformational space in solution.

Open image in new window

Graphical Abstract Emergence of the Bcl-2-rFLD’s structural heterogeneity in solution, evidenced by molecular dynamics simulation.

     

Multiple loop conformations of peptides predicted by molecular dynamics simulations are compatible with nuclear magnetic resonance

The affinity and selectivity of protein-protein interactions can be fine-tuned by varying the size, flexibility, and amino acid composition of involved surface loops. As a model for such surface loops, we study the conformational landscape of an octapeptide, whose flexibility is chemically steered by a covalent ring closure integrating an azobenzene dye into and by a disulfide bridge additionally constraining the peptide backbone. Because the covalently integrated azobenzene dyes can be switched by light between a bent cis state and an elongated trans state, six cyclic peptide models of strongly different flexibilities are obtained. The conformational states of these peptide models are sampled by NMR and by unconstrained molecular dynamics (MD) simulations. Prototypical conformations and the free-energy landscapes in the high-dimensional space spanned by the phi/psi angles at the peptide backbone are obtained by clustering techniques from the MD trajectories. Multiple open-loop conformations are shown to be predicted by MD particularly in the very flexible cases and are shown to comply with the NMR data despite the fact that such open-loop conformations are missing in the refined NMR structures.     

Structural elucidation of 4'-epiadriamycin by nuclear magnetic resonance spectroscopy and comparison with adriamycin and daunomycin using quantum mechanical and restrained molecular dynamics approach

The structural and electronic properties of 4′-epiadriamycin, adriamycin, and daunomycin have been studied using density functional theory (DFT) employing B3LYP exchange correlation. The chemical shifts of ¹H and ¹³C resonances in nuclear magnetic resonance spectra have been calculated using Gauge-Invariant Atomic Orbital (GIAO) method as implemented in Gaussian 98 and compared with experimental spectra recorded at 500 MHz. ¹³C resonances of drugs have been assigned for the first time. A restrained molecular dynamics approach was used to get the optimized solution structure of drugs using inter-proton distance constraints obtained from 2D NOESY spectra. The glycosidic angle C7-O7-C1′-C2′ is found to show considerable flexibility by adopting 156°-161° (I), 142°-143° (II), and 38°-78° (III) conformations, of which the biological relevant structure appears to be the conformer II. The observed different conformations of the three drugs are correlated to the differential anticancer activity and the available biochemical evidence exhibited by these drugs.     

Hairpin structures in DNA containing arabinofuranosylcytosine. A combination of nuclear magnetic resonance and molecular dynamics

Nuclear magnetic resonance (NMR) and model-building studies were carried out on the hairpin form of the octamer d(CGaCTAGCG) (aC = arabinofuranosylcytosine), referred to as the TA compound. The nonexchangeable protons of the TA compound were assigned by means of nuclear Overhauser effect spectroscopy (NOESY) and correlated spectroscopy (COSY). From a detailed analysis of the coupling data and of the NOESY spectra the following conclusions are reached: (i) The hairpin consists of a stem of three Watson-Crick type base pairs, and the two remaining residues, T(4) and dA(5), participate in a loop. (ii) All sugar rings show conformational flexibility although a strong preference for the S-type (C2'-endo) conformer is observed. (iii) The thymine does not stack upon the 3' side of the stem as expected, but swings into the minor groove. (This folding principle of the loop involves an unusual alpha t conformer in residue T(4).) (iv) At the 5'-3' loop-stem junction a stacking discontinuity occurs as a consequence of a sharp turn in that part of the backbone, caused by the unusual beta + and gamma t torsion angles in residue dG(6). (v) The A base slides over the 5' side of the stem to stack upon the aC(3) residue at the 3' side of the stem in an antiparallel fashion. On the basis of J couplings and a set of approximate proton-proton distances from NOE cross peaks, a model for the hairpin was constructed. This model was then refined by using an iterative relaxation matrix approach (IRMA) in combination with restrained molecular dynamics calculations. The resulting final model satisfactorily explains all the distance constraints.     

Structural characterization of cyclic kallidin analogues in DMSO by nuclear magnetic resonance and molecular dynamics.

The conformational properties in DMSO of two head-to-tail cyclic analogues of kallidin ([Lys(0)]-bradykinin, KL) as well as those of the corresponding linear peptides were studied by NMR and molecular dynamics (MD) simulations. The modifications in the sequence were introduced at position 6, resulting in the four peptides, [Tyr(6)]-KL (YKL), [Trp(6)]-KL (WKL), cyclo-([Tyr(6)]-KL) (YCKL) and cyclo-([Trp(6)]-KL) (WCKL).The linear WKL analogue was significantly more potent than kallidin on rat duodenum preparations, whereas YKL was significantly less potent. Both cyclic peptides, YCKL and WCKL displayed similar activity, lower than that of the linear analogues and also of cyclo-KL.The two linear analogues display high conformational flexibility in DMSO. In the predominant conformer, for both peptides, all three X-Pro bonds adopt a trans configuration. Three out of four conformers present in YCKL and WCKL were completely assigned. The configurations at the X-Pro bonds are the same for the two analogues. All cyclic conformers show a cis configuration in at least one X-Pro bond and always opposite configuration for the two consecutive X-Pro bonds.The NOE-restrained MD calculations resulted in the detection of several elements of secondary structure in each of the conformers. Such elements are described and their possible relevance to biological activity is discussed.     

Structure of DNA sequence d-TGATCA by two-dimensional nuclear magnetic resonance spectroscopy and restrained molecular dynamics

The 5' d-TpG 3' element is a part of DNA sequences involved in regulation of gene expression and is also a site for intercalation of several anticancer drugs. Solution conformation of DNA duplex d-TGATCA containing this element has been investigated by two-dimensional NMR spectroscopy. Using a total of 12 torsional angles and 121 distance constraints, structural refinement has been carried out by restrained molecular dynamics (rMDs) in vacuum up to 100 ps. The structure is characterized by a large positive roll at TpG/CpA base pair step and large negative propeller twist for AT and TA base pairs. The backbone torsional angle, gamma(O5'-C5'-C4'-C3'), of T1 residue adopts a trans-conformation which is corroborated by short intra nucleotide T1H6-T1H5' (3.7A) distance in nuclear overhauser effect spectroscopy (NOESY) spectra while the backbone torsional angle, beta(P-O5'-C5'-C4'), exists in trans as well as gauche state for T1 and C5 residues. There is evidence of significant flexibility of the sugar-phosphate backbone with rapid inter-conversion between two different conformers at TpG/CpA base pair step. The base sequence dependent variations and local structural heterogeneity have important implications in specific recognition of DNA by ligands.     

Molecular dynamics simulation of hen egg white lysozyme: a test of the GROMOS96 force field against nuclear magnetic resonance data

Biomolecular force fields for use in molecular dynamics (MD) simulations of proteins, DNA, or membranes are generally parametrized against ab initio quantum-chemical and experimental data for small molecules. The application of a force field in a simulation of a biomolecular system, such as a protein in solution, may then serve as a test of the quality and transferability of the force field. Here, we compare various properties obtained from two MD simulations of the protein hen egg white lysozyme (HEWL) in aqueous solution using the latest version, GROMOS96, of the GROMOS force field and an earlier version, GROMOS87+, with data derived from nuclear magnetic resonance (NMR) experiments: NOE atom-atom distance bounds, (3)J(HNalpha)-coupling constants, and backbone and side-chain order parameters. The convergence of these quantities over a 2-ns period is considered, and converged values are compared to experimental ones. The GROMOS96 simulation shows better agreement with the NMR data and also with the X-ray crystal structure of HEWL than the GROMOS87+ simulation, which was based on an earlier version of the GROMOS force field.     

Insights into internal dynamics of 6-phosphogluconolactonase from Trypanosoma brucei studied by nuclear magnetic resonance and molecular dynamics

Nuclear magnetic resonance is used to investigate the backbone dynamics in 6-phosphogluconolactonase from Trypanosoma brucei (Tb6PGL) with (holo-) and without (apo-) 6-phosphogluconic acid as ligand. Relaxation data were analyzed using the model-free approach and reduced spectral density mapping. Comparison of predictions, based on 77 ns molecular dynamics simulations, with the observed relaxation rates gives insight into dynamical properties of the protein and their alteration on ligand binding. Data indicate dynamics changes in the vicinity of the binding site. More interesting is the presence of perturbations located in remote regions of this well-structured globular protein in which no large-amplitude motions are involved. This suggests that delocalized changes in dynamics that occur upon binding could be a general feature of protein-target interactions.     

Protein structure and dynamics in nonaqueous solvents: insights from molecular dynamics simulation studies

Protein structure and dynamics in nonaqueous solvents are here investigated using molecular dynamics simulation studies, by considering two model proteins (ubiquitin and cutinase) in hexane, under varying hydration conditions. Ionization of the protein groups is treated assuming "pH memory," i.e., using the ionization states characteristic of aqueous solution. Neutralization of charged groups by counterions is done by considering a counterion for each charged group that cannot be made neutral by establishing a salt bridge with another charged group; this treatment is more physically reasonable for the nonaqueous situation, contrasting with the usual procedures. Our studies show that hydration has a profound effect on protein stability and flexibility in nonaqueous solvents. The structure becomes more nativelike with increasing values of hydration, up to a certain point, when further increases render it unstable and unfolding starts to occur. There is an optimal amount of water, approximately 10% (w/w), where the protein structure and flexibility are closer to the ones found in aqueous solution. This behavior can explain the experimentally known bell-shaped dependence of enzyme catalysis on hydration, and the molecular reasons for it are examined here. Water and counterions play a fundamental and dynamic role on protein stabilization, but they also seem to be important for protein unfolding at high percentages of bound water.     

Thermal unfolding molecular dynamics simulation of Escherichia coli dihydrofolate reductase: thermal stability of protein domains and unfolding pathway

Temperature induced unfolding of Escherichia coli dihydrofolate reductase was carried out by using molecular dynamic simulations. The simulations show that the unfolding generally involves an initial end-to-end collapse of the adenine binding domain into partially extended loops, followed by a gradual breakdown of the remaining beta sheet core structure. The core, which consists of beta strands 5-7, was observed to be the most resistant to thermal unfolding. This region, which is made up of part of the N terminus domain and part of the large domain of the E. coli dihydrofolate reductase, may constitute the nucleation site for protein folding and may be important for the eventual formation of both domains. The unfolding of different domains at different stages of the unfolding process suggests that protein domains vary in stability and that the rate at which they unfold can affect the overall outcome of the unfolding pathway. This observation is compared with the recently proposed hierarchical folding model. Finally, the results of the simulation were found to be consistent with a previous experimental study (Frieden, Proc Natl Acad Sci USA 1990;87:4413-4416) which showed that the folding process of E. coli dihydrofolate reductase involves sequential formation of the substrate binding sites.     

Protein dynamics by solid-state nuclear magnetic resonance spectroscopy: Peptide backbone of the coat protein in fd bacteriophage

¹³C and ¹⁵N chemical shift anisotropy and ¹⁵N¹H dipolar powder patterns from backbone sites of the coat protein in fd bacteriophage are not averaged by motion. This means that the polypeptide backbone of the protein has no large amplitude motions rapid compared to 10⁴ Hz. Relaxation studies on the ¹³C_α and ¹⁵N amide resonances indicate the presence of motions on the 10⁹ Hz timescale. These results are reconciled with a model where an otherwise rigid backbone undergoes small amplitude, rapid motions.     

A comparison of 15N NMR relaxation measurements with a molecular dynamics simulation: Backbone dynamics of the glucocorticoid receptor DNA-binding domain

The rapid motions of the backbone of the DNA-binding domain of the glucocorticoid receptor (GR DBD) have been investigated using proton-detected heteronuclear NMR experiments on ¹⁵N-labeled protein at pH 6.0 and with a 200 psec molecular dynamics simulation of hydrated GR DBD. The experimental data were interpreted in terms of a generalized order parameter (S²) and an effective correlation time (τ_e) for the internal motion of each amide bond. A back calculation, using the same model, yielded the {¹H}-¹⁵N nuclear Overhauser effects (NOEs) and the ¹⁵N spin-lattice relaxation times (T₁) from the simulated data. The rapid motions of the backbone turned out to be rather limited and uniform throughout the protein, with a somewhat reduced mobility in the two major α-helical regions and a slightly enhanced flexibility for some residues in the first zinc coordinating region. The agreement between the experimental and simulated S²-values was as good as quantitative for most of the residues, except for some residues that were subject to a more large-scale, and in the simulation thus poorly sampled, motion. Examples of such motions that were found in the simulation include jumps of the amide bond of Ile-487 between the charged oxygens of the side chain of Asp-485 and less distinct large scale motions for some of the residues in the extended regions, that were shown to give rise to noisy and/or fast decaying internal reorientational correlation functions. For these residues large differences in the simulated and experimental τ_e-values were found, indicating that motions on different time scales were dominating in the experimental and simulated values. The lower (<0.7) experimental NOEs for these residues could not be reproduced in the simulation and were shown to be a consequence of the lower τ_e-values estimated in the simulation. By combining information from the simulation and the experiment a more complete picture of the motions for these residues can be obtained as is illustrated with an estimation of the jump angle and jump frequency for the amide bond of Ile-487. © 1993 Wiley-Liss, Inc.