Pathways to a protein folding intermediate observed in a 1-microsecond simulation in aqueous solution.

An implementation of classical molecular dynamics on parallel computers of increased efficiency has enabled a simulation of protein folding with explicit representation of water for 1 microsecond, about two orders of magnitude longer than the longest simulation of a protein in water reported to date. Starting with an unfolded state of villin headpiece subdomain, hydrophobic collapse and helix formation occur in an initial phase, followed by conformational readjustments. A marginally stable state, which has a lifetime of about 150 nanoseconds, a favorable solvation free energy, and shows significant resemblance to the native structure, is observed; two pathways to this state have been found.

Calculation of the relative change in binding free energy of a protein-inhibitor complex.

By means of a thermodynamic perturbation method implemented with molecular dynamics, the relative free energy of binding was calculated for the enzyme thermolysin complexed with a pair of phosphonamidate and phosphonate ester inhibitors. The calculated difference in free energy of binding was 4.21 +/- 0.54 kilocalories per mole. This compares well with the experimental value of 4.1 kilocalories per mole. The method is general and can be used to determine a change or "mutation" in any system that can be suitably represented. It is likely to prove useful for protein and drug design.

Free energy calculations by computer simulation.

A fundamental problem in chemistry and biochemistry is understanding the role of solvation in determining molecular properties. Recent advances in statistical mechanical theory and molecular dynamics methodology can be used to solve this problem with the aid of supercomputers. By using these advances the free energies of solvation of all the chemical classes of amino acid side chains, four nucleic acid bases and other organic molecules can be calculated. The effect of a site-specific mutation on the stability of trypsin is predicted. The results are in good agreement with available experiments.

Probing protein stability with unnatural amino acids.

Unnatural amino acid mutagenesis, in combination with molecular modeling and simulation techniques, was used to probe the effect of side chain structure on protein stability. Specific replacements at position 133 in T4 lysozyme included (i) leucine (wt), norvaline, ethylglycine, and alanine to measure the cost of stepwise removal of methyl groups from the hydrophobic core, (ii) norvaline and O-methyl serine to evaluate the effects of side chain solvation, and (iii) leucine, S,S-2-amino-4-methylhexanoic acid, and S-2-amino-3-cyclopentylpropanoic acid to measure the influence of packing density and side chain conformational entropy on protein stability. All of these factors (hydrophobicity, packing, conformational entropy, and cavity formation) significantly influence protein stability and must be considered when analyzing any structural change to proteins.

Molecular dynamics studies of the HIV-1 TAR and its complex with argininamide.

The dynamic behavior of HIV-1 TAR and its complex with argininamide is investigated by means of molecular dynamics simulations starting from NMR structures, with explicit inclusion of water and periodic boundary conditions particle mesh Ewald representation of the electrostatic energy. During simulations of free and argininamide-bound TAR, local structural patterns, as determined by NMR experiments, were reproduced. An interdomain motion was observed in the simulations of free TAR, which is absent in the case of bound TAR, leading to the conclusion that the free conformation of TAR is intrinsically more flexible than the bound conformation. In particular, in the bound conformation the TAR-argininamide interface is very well ordered, as a result of the formation of a U.A.U base triple, which imposes structural constraints on the global conformation of the molecule. Free energy analysis, which includes solvation contributions, was used to evaluate the influence of van der Waals and electrostatic terms on formation of the complex and on the conformational rearrangement from free to bound TAR.

Insight into the stabilization of A-DNA by specific ion association: spontaneous B-DNA to A-DNA transitions observed in molecular dynamics simulations of d[ACCCGCGGGT]2 in the presence of hexaamminecobalt(III).

BACKGROUND: Duplex DNA is more than a simple information carrier. The sequence-dependent structure and its inherent deformability, in concert with the subtle modulating effects of the environment, play a crucial role in the regulation and packaging of DNA. Recent advances in force field and simulation methodologies allow molecular dynamics simulations to now represent the specific effects of the environment. An understanding of the environmental dependence of DNA structure gives insight into how histones are able to package DNA, how various proteins are able to bind and modulate nucleic acid structure and will ultimately aid the design of molecules to package DNA for more effective gene therapy. RESULTS: Molecular dynamics simulations of d[ACCCGCGGGT]2 in solution in the presence of hexaamminecobalt(III) [Co(NH3)6(3+)] show stabilization of A-DNA and spontaneous B-DNA to A-DNA transitions, which is consistent with experimental results from NMR and Raman spectroscopic and X-ray crystallographic studies. In the absence of Co(NH3)6(3+), A-DNA to B-DNA transitions are observed instead. In addition to their interaction with the guanines in the major groove, Co(NH3)6(3+) ions bridge opposing strands in the bend across the major groove, probably stabilizing A-DNA. CONCLUSIONS: The simulation methods and force fields have advanced to a sufficient level that some representation of the environment can be seen in nanosecond length molecular dynamics simulations. These simulations suggest that, in addition to the general explanation of A-DNA stabilization by dehydration, hydration and ion association in the major groove stabilize A-DNA.

Free-energy calculations highlight differences in accuracy between X-ray and NMR structures and add value to protein structure prediction.

BACKGROUND: While X-ray crystallography structures of proteins are considerably more reliable than those from NMR spectroscopy, it has been difficult to assess the inherent accuracy of NMR structures, particularly the side chains. RESULTS: For 15 small single-domain proteins, we used a molecular mechanics-/dynamics-based free-energy approach to investigate native, decoy, and fully extended alpha conformations. Decoys were all less energetically favorable than native conformations in nine of the ten X-ray structures and in none of the five NMR structures, but short 150 ps molecular dynamics simulations on the experimental structures caused them to have the lowest predicted free energy in all 15 proteins. In addition, a strong correlation exists (r(2) = 0.86) between the predicted free energy of unfolding, from native to fully extended conformations, and the number of residues. CONCLUSIONS: This work suggests that the approximate treatment of solvent used in solving NMR structures can lead NMR model conformations to be less reliable than crystal structures. This conclusion was reached because of the considerably higher calculated free energies and the extent of structural deviation during aqueous dynamics simulations of NMR models compared to those determined by X-ray crystallography. Also, the strong correlation found between protein length and predicted free energy of unfolding in this work suggests, for the first time, that a free-energy function can allow for identification of the native state based on calculations on an extended state and in the absence of an experimental structure.

Molecular dynamics study of the stability of staphylococcal nuclease mutants: component analysis of the free energy difference of denaturation.

The stability of two mutants G88V (Gly-88-->Val) and A69T (Ala-69-->Thr) of staphylococcal nuclease was analyzed by molecular dynamics simulations. The calculated free energy differences of denaturation for G88V and A69T were -1.1 and -2.8 kcal/mol, respectively. These values are in good agreement with the experimental values. The free energy differences divided into electrostatic and van der Waals components were analyzed. These two mutants are mainly destabilized due to van der Waals interactions. There is little difference between the electrostatic contribution to the free energy change in the native state and that in the denatured state. In each mutant structure, a small cavity appears in the vicinity of the perturbed residue. It is suggested that intramolecular van der Waals interactions of the mutants are weaker than those of the wild-type. Furthermore, analyses of the contributions of each residue near the perturbed residue and of water to the free energy difference of denaturation suggest that the interaction between water and the perturbed residue plays a very important role in the stability of staphylococcal nuclease, and that a small hydrophobic core consisting of the three aromatic rings (Tyr-27, Phe-34, Phe-76) and the side chain of Met-32 is also important for the stability.

A conformational analysis of daunomycin considering D-ring repucker.

We present a complete molecule mechanics optimization of daunomycin. Two different D-ring puckers are found to be of comparable energy, consistent with NMR data, although only one of these low-energy structures has been observed by X-ray crystallography. Our results are more consistent with the previous conformational analysis of daunomycin by Neidle and Taylor (Neidle, S. and Taylor, G.L. (1979) FEBS Lett. 107, 348-354) than that of Nakata and Hopfinger (Nakata, Y. and Hopfinger, A.J. (1980) FEBS Lett. 117, 259-264).

An NMR study of the interaction of daunomycin with dinucleotides and dinucleoside phosphates.

The mode of action of the widely used anticancer drug daunomycin was studied by 360 MHz and 100 MHz proton NMR. Information obtained from scalar coupling constants indicates that the conformation of the A ring of the aglycone moiety of daunomycin is maintained upon binding to dinucleotides but that the conformation of the daunosamine ring moiety is altered upon binding. Ring-current-induced chemical shifts of the drug protons were observed as the drug was titrated with deoxydinucleotides. The chemical shift results that daunomycin forms 1 : 1 complexes with deoxydinucleotides and that the strength of the drug-nucleotide interaction is dependent upon the base composition of the dinucleotide. On the basis of the NMR data, a model for the drug/deoxydinucleotide complex is presented featuring intercalation of the drug between the bases of the dinucleoside.

Investigating the binding specificity of U1A-RNA by computational mutagenesis.

The mammalian spliceosomal protein U1A binds a hairpin RNA with picomolar affinity. To examine the origin of this binding specificity, we carried out computational mutagenesis on protein and RNA residues in the U1A-RNA binding interface. Our computational mutagenesis methods calculate the relative binding affinity between mutant and wild-type as the sum of molecular mechanical energies and solvation free energies estimated with a continuum solvent model. We obtained good agreement with experimental studies and we verified mutations that abolish and improve binding. Therefore, we offer these methods as computationally inexpensive tools for investigating and predicting the effects of site-specific mutagenesis.

Calculation of the relative binding free energy of 2'GMP and 2'AMP to ribonuclease T1 using molecular dynamics/free energy perturbation approaches.

We present a calculation of the relative changes in binding free energy between the complex of ribonuclease T1 (RNase Tr) with its inhibitor 2'-guanosine monophosphate (2'GMP) and that of RNase T1-2'-adenosine monophosphate (2'AMP) by means of a thermodynamic perturbation method implemented with molecular dynamics. Using the available crystal structure of the RNase T1-2'GMP complex, the structure of the RNase T1-2'AMP complex was obtained as a final structure of the perturbation calculation. The calculated difference in the free energy of binding (delta delta Gbind) was 2.76 kcal/mol. This compares well with the experimental value of 3.07 kcal/mol. The encouraging agreement in delta delta Gbind suggests that the interactions of inhibitors with the enzyme are reasonably represented. Energy component analyses of the two complexes reveal that the active site of RNase T1 electrostatically stabilizes the binding of 2'GMP more than that of 2'AMP by 44 kcal/mol, while the van der Waals' interactions are similar in the two complexes. The analyses suggest that the mutation from Glu46 to Gln may lead to a preference of RNase T1 for adenine in contrast to the guanine preference of the wild-type enzyme. Although the molecular dynamics equilibration moves the atoms of the RNase T1-2'GMP system about 0.9 A from their X-ray positions and the mutation of the G to A in the active site increases the deviation from the X-ray structure, the mutation of the A back to G reduces the deviation. This and the agreement found for delta delta Gbind suggest that the molecular dynamics/free energy perturbation method will be useful for both energetic and structural analysis of protein-ligand interactions.

Free energy calculations on dimer stability of the HIV protease using molecular dynamics and a continuum solvent model.

Dimerization of HIV-I protease (HIV PR) monomers is an essential prerequisite for viral proteolytic activity and the subsequent generation of infectious virus particles. Disrupting dimerization of the enzyme can inhibit its activity. We have calculated the relative binding free energies between different dimers of the HIV protease using molecular dynamics and a continuum model, which we call MM/PBSA. We examined the dominant negative inhibition of the HIV PR by a mutated form of the protease and found relative dimerization free energies of homo- and hetero-dimerization consistent with experimental data. We also developed a rapid screening method, which was called the virtual mutagenesis method to consider other mutations which might stabilize non-wild-type heterodimers. Using this approach, we considered the mutations near the dimer interface which might cause dominant negative inhibition of the HIV PR. The rapid method we developed can be used in studying any ligand-protein and protein-protein interaction, in order to identify mutations that can enhance the binding affinities of the complex.

Are time-averaged restraints necessary for nuclear magnetic resonance refinement? A model study for DNA.

A recently suggested method for refinement of structural data obtained from two-dimensional nuclear magnetic resonance experiments using molecular dynamics (MD) is explored. In this method, the time-averaged values of the appropriate internal co-ordinates of the molecule, calculated from the MD trajectory, are driven by restraints towards the experimental target values. This contrasts with most refinement procedures currently in use, where restraints are applied based on the instantaneous values of the appropriate co-ordinates. Both refinement methods are applied to the EcoRI restriction site DNA hexamer d(GAATTC)2, using target nuclear Overhauser enhancement distances derived from a one nanosecond unrestrained MD simulation of this structure. The resulting refined structures are compared to the results of the unrestrained MD trajectory, which serves as our "experimental" data. We show that although both methods can yield an average structure with the correct gross morphology, the new method allows both a much more realistic picture of inherent flexibility, and reproduces fine conformational detail better, such as sequence dependency. We also analyze the very long MD trajectory generated here (longer than any previously reported for a DNA oligomer), and find that significantly shorter simulations, typical of those frequently performed, may not yield acceptably reliable values for certain structural parameters.

Theoretical studies of an exceptionally stable RNA tetraloop: observation of convergence from an incorrect NMR structure to the correct one using unrestrained molecular dynamics.

We report on the results of five independent and unrestrained molecular dynamics simulations of an RNA tetraloop, r(GGACUUCGGUCC), and its related structures with the loop UUCG sugars changed to deoxyribose. Two separate NMR structures have been reported for the loop portion of this molecule, with the second refinement resulting in a slightly different and more accurate conformation for the loop. The root-mean-square deviation (RMSd) between the two NMR structures, for the loop portions only, is 2.5 A. Our simulations, starting from the two NMR structures, demonstrate that this tetraloop is a very stable and rigid structure with both nanosecond length simulations staying very close to the initial structures. Additionally, both simulations preserved most, if not all, of the NMR-derived interactions and violated very few of the nuclear Overhauser effect (NOE)-derived distances used in the structure refinements. However, when the two NMR structures were simulated with deoxyriboses in the loops instead of the native riboses, the flexibility of the systems increased and we observed a conversion from the incorrect to the correct loop conformation in the simulation which started in the incorrect loop conformation. When the riboses were subsequently re-introduced back into the structure which underwent the conversion, the agreement between this simulation and the one starting from the correct NMR structure was a remarkably low 0.5 A, demonstrating an almost complete convergence from the incorrect to the correct structure using unrestrained molecular dynamics.

Structure and thermodynamics of RNA-protein binding: using molecular dynamics and free energy analyses to calculate the free energies of binding and conformational change.

An adaptive binding mechanism, requiring large conformational rearrangements, occurs commonly with many RNA-protein associations. To explore this process of reorganization, we have investigated the conformational change upon spliceosomal U1A-RNA binding with molecular dynamics (MD) simulations and free energy analyses. We computed the energetic cost of conformational change in U1A-hairpin and U1A-internal loop binding using a hybrid of molecular mechanics and continuum solvent methods. Encouragingly, in all four free energy comparisons (two slightly different proteins, two different RNAs), the free macromolecule was more stable than the bound form by the physically reasonable value of approximately 10 kcal/mol. We calculated the absolute binding free energies for both complexes to be in the same range as that found experimentally.

Observation of the A-DNA to B-DNA transition during unrestrained molecular dynamics in aqueous solution.

A large challenge in molecular dynamics (MD) simulations of proteins and nucleic acids is to find the correct "experimental" geometry when a simulation is started a significant distance away from it. In this study, we have carried out four unrestrained approximately 1 ns length MD trajectories in aqueous solution on the DNA duplex d(CCAACGTTGG)2, two beginning in a canonical A-DNA structure and two beginning in a canonical B-DNA structure. As judged by root-mean-squared coordinate deviations, average structures computed from all four of the trajectories converge to within approximately 0.8 to 1.6 angstroms (all atoms) of each other, which is 1.3 to 1.7 angstroms (all atoms of the central six residues from each strand) and 3.1 to 3.6 angstroms (all atoms) away from the B-DNA-like X-ray structure reported for this sequence. To our knowledge, this is the first example of multiple nanosecond molecular dynamics trajectories with full representation of DNA charges, solvent and long range electrostatics that demonstrate both internal consistency (two different starting structures and four different trajectories lead to a consistent average structure) and considerable agreement with the X-ray crystal structure of this sequence and NMR data on duplex DNA in aqueous solution. This internal consistency of structure for a given sequence suggests that one can now begin to realistically examine sequence-dependent structural effects in DNA duplexes using molecular dynamics.

Molecular dynamics simulations of "loop closing" in the enzyme triose phosphate isomerase.

We present molecular dynamics simulations on the active site region of dimeric triose phosphate isomerase (TIM) using the co-ordinates of native chicken muscle TIM as a starting point and performing simulations with no substrate, with dihydroxyacetone phosphate (DHAP), the natural substrate, and with dihydroxyacetone sulfate (DHAS), a substrate analog. Whereas most of the protein moves less than 1 A during the simulation, some residues in the active site loop move more than 8 A during the 10.5 picoseconds of dynamics for each of the simulations. Most interestingly, the nature of the loop motion depends on the substrate, with the largest motion found in the presence of DHAP, and only in the presence of DHAP does the loop move to "close off" the active site pocket. The final structure found for the DHAP-chicken TIM complex is qualitatively similar to that described by Alber et al. for DHAP-yeast TIM. Simulations on the monomeric protein gives insight into why the molecule is active only as a dimer.

Simple model for the effect of Glu165----Asp165 mutation on the rate of catalysis in triose phosphate isomerase.

We present an ab-initio self-consistent field calculation with a 4-31G basis set on a simple model for proton abstraction from hydroxyacetone (a model for dihydroxyacetone phosphate; DHAP) by formate, which is a model for Glu165 in triose phosphate isomerase. Earlier, we showed that the electrophilic groups on the enzyme (the NH3+ of Lys13 and the NH of His95) were essential to efficient catalysis by triose phosphate isomerase. These groups stabilized the enediolate formed by proton abstraction from the DHAP model so that proton transfer from this molecule to Glu165 became likely. In this study, we carry this analysis one step further. First, we re-examine the energy profile for proton transfer, using the fact that our earlier calculations showed that the combined effect of His95 and Lys13 on the reactant DHAP and intermediate enediolate was to make them equal in energy. Then, we analyze the likely effect of changing Glu165 to Asp165 and relate this to experiments on the kinetics of enzyme catalysis by the Glu165----Asp165 mutant.

Molecular dynamics in the endgame of protein structure prediction.

In order adequately to sample conformational space, methods for protein structure prediction make necessary simplifications that also prevent them from being as accurate as desired. Thus, the idea of feeding them, hierarchically, into a more accurate method that samples less effectively was introduced a decade ago but has not met with more than limited success in a few isolated instances. Ideally, the final stages should be able to identify the native state, show a good correlation with native similarity in order to add value to the selection process, and refine the structures even further. In this work, we explore the possibility of using state-of-the-art explicit solvent molecular dynamics and implicit solvent free energy calculations to accomplish all three of those objectives on 12 small, single-domain proteins, four each of alpha, beta and mixed topologies. We find that this approach is very successful in ranking the native and also enhances the structure selection of predictions generated from the Rosetta method.