Protein conformational dynamics of crambin in crystal, solution and in the trajectories of molecular dynamics simulations

Atomic displacement parameters — B factors of the eight crambin crystal structures obtained at 0.54–1.5 Å resolution and temperatures of 100–293 K have been analyzed. The comparable contributions to the B factor values are the intramolecular motions which are modeled by the harmonic vibration calculations and derived from the molecular dynamics simulation (MD) as well as rigid body changes in the position of a protein molecule as a whole. In solution for the average NMR structure of crambin the amplitudes of the backbone atomic fluctuations of the most residues of the segments with the regular backbone conformations are close to the amplitudes of the small scale harmonic vibrations. For the same residues the probability of the medium scale fluctuations fixed by the hydrogen exchange method is very low. The restricted conformational mobility of those segments is coupled with the depressed amplitudes of the fluctuation changes of the tertiary structure registered by the residue accessibility changes in an ensemble of NMR structures that forms the average NMR structure of crambin. The amplitudes of temperature fluctuations of backbone atoms and the tertiary structure raise in the segment with the irregular conformations, turn and loops. In the same segments the amplitudes of the calculated harmonic vibrations also increase, but to a lesser extent and especially in the interhelical loop with the most strong and complicated fluctuation changes of the backbone conformation. In solution for the NMR structure in this loop the conformational transitions occur between the conformational substates separated by the energy barriers, but they are not observed even in the long 100 ns trajectories from the MD simulation of crambin. These strong local fluctuation changes of the structure may play a key role in the protein functioning and modern performance improvements in the MD simulation techniques are oriented to increase the probability of protein appearance in the trajectories from the MD simulations.     

Locally accessible conformations of proteins: multiple molecular dynamics simulations of crambin.

Multiple molecular dynamics (MD) simulations of crambin with different initial atomic velocities are used to sample conformations in the vicinity of the native structure. Individual trajectories of length up to 5 ns sample only a fraction of the conformational distribution generated by ten independent 120 ps trajectories at 300 K. The backbone atom conformational space distribution is analyzed using principal components analysis (PCA). Four different major conformational regions are found. In general, a trajectory samples only one region and few transitions between the regions are observed. Consequently, the averages of structural and dynamic properties over the ten trajectories differ significantly from those obtained from individual trajectories. The nature of the conformational sampling has important consequences for the utilization of MD simulations for a wide range of problems, such as comparisons with X-ray or NMR data. The overall average structure is significantly closer to the X-ray structure than any of the individual trajectory average structures. The high frequency (less than 10 ps) atomic fluctuations from the ten trajectories tend to be similar, but the lower frequency (100 ps) motions are different. To improve conformational sampling in molecular dynamics simulations of proteins, as in nucleic acids, multiple trajectories with different initial conditions should be used rather than a single long trajectory.     

Structure and internal dynamics of the bovine pancreatic trypsin inhibitor in aqueous solution from long-time molecular dynamics simulations

Structural and dynamic properties of bovine pancreatic trypsin inhibitor (BPTI) in aqueous solution are investigated using two molecular dynamics (MD) simulations: one of 1.4 ns length and one of 0.8 ns length in which atom-atom distance bounds derived from NMR spectroscopy are included in the potential energy function to make the trajectory satisfy these experimental data more closely. The simulated properties of BPTI are compared with crystal and solution structures of BPTI, and found to be in agreement with the available experimental data. The best agreement with experiment was obtained when atom-atom distance restraints were applied in a time-averaged manner in the simulation. The polypeptide segments found to be most flexible in the MD simulations coincide closely with those showing differences between the crystal and solution structures of BPTI. © 1995 Wiley-Liss, Inc.     

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.     

Using molecular dynamics simulations on crambin to evaluate the suitability of different continuum dielectric and hydrogen atom models for protein simulations

Molecular dynamics simulations of enzymes with enough explicit waters of solvation to realistically account for solute-solvent interactions can burden the computational resources required to perform the simulation by more than two orders of magnitude. Since enzyme simulations even with an implicit solvation model can be imposing for a supercomputer, it is important to assess the suitability of different continuum dielectric models for protein simulations. A series of 100-picosecond molecular dynamics simulations were performed on the X-ray crystal structure of the protein crambin to examine how well computed structures, obtained using seven continuum dielectric and two hydrogen atom models, agreed with the X-ray structure. The best level of agreement between computed and experimental structures was obtained using a constant dielectric of 2 and the all-hydrogen model. Continuum dielectric models of 1, 1r, and 2r also led to computed structures in reasonably good agreement with the X-ray structure. In all cases, the all-hydrogen model gave better agreement than the united atom model, although, in one case, the difference was not significant. Dielectric models of 4, 80, and 4r with either hydrogen model yielded significantly poorer fits. It is especially noteworthy that the observed trends did not semiquantitatively converge until about 50 picoseconds into the simulations, suggesting that validation studies for protein calculations based on energy minimizations or short simulations should be viewed with caution.     

Structure, dynamics and hydration of the nogalamycin-d(ATGCAT)2Complex determined by NMR and molecular dynamics simulations in solution.

The structure of the 1:1 nogalamycin:d(ATGCAT)2 complex has been determined in solution from high-resolution NMR data and restrained molecular dynamics (rMD) simulations using an explicit solvation model. The antibiotic intercalates at the 5'-TpG step with the nogalose lying along the minor groove towards the centre of the duplex. Many drug-DNA nuclear Overhauser enhancements (NOEs) in the minor groove are indicative of hydrophobic interactions over the TGCA sequence. Steric occlusion prevents a second nogalamycin molecule from binding at the symmetry-related 5'-CpA site, leading to the conclusion that the observed binding orientation in this complex is the preferred orientation free of the complication of end-effects (drug molecules occupy terminal intercalation sites in all X-ray structures) or steric interactions between drug molecules (other NMR structures have two drug molecules bound in close proximity), as previously suggested. Fluctuations in key structural parameters such as rise, helical twist, slide, shift, buckle and sugar pucker have been examined from an analysis of the final 500 ps of a 1 ns rMD simulation, and reveal that many sequence-dependent structural features previously identified by comparison of different X-ray structures lie within the range of dynamic fluctuations observed in the MD simulations. Water density calculations on MD simulation data reveal a time-averaged pattern of hydration in both the major and minor groove, in good agreement with the extensive hydration observed in two related X-ray structures in which nogalamycin is bound at terminal 5'-TpG sites. However, the pattern of hydration determined from the sign and magnitude of NOE and ROE cross-peaks to water identified in 2D NOESY and ROESY experiments identifies only a few "bound" water molecules with long residence times. These solvate the charged bicycloaminoglucose sugar ring, suggesting an important role for water molecules in mediating drug-DNA electrostatic interactions within the major groove. The high density of water molecules found in the minor groove in X-ray structures and MD simulations is found to be associated with only weakly bound solvent in solution.     

Simulation of nitroxide electron paramagnetic resonance spectra from brownian trajectories and molecular dynamics simulations

A simulated continuous wave electron paramagnetic resonance spectrum of a nitroxide spin label can be obtained from the Fourier transform of a free induction decay. It has been previously shown that the free induction decay can be calculated by solving the time-dependent stochastic Liouville equation for a set of Brownian trajectories defining the rotational dynamics of the label. In this work, a quaternion-based Monte Carlo algorithm has been developed to generate Brownian trajectories describing the global rotational diffusion of a spin-labeled protein. Also, molecular dynamics simulations of two spin-labeled mutants of T4 lysozyme, T4L F153R1, and T4L K65R1 have been used to generate trajectories describing the internal dynamics of the protein and the local dynamics of the spin-label side chain. Trajectories from the molecular dynamics simulations combined with trajectories describing the global rotational diffusion of the protein are used to account for all of the dynamics of a spin-labeled protein. Spectra calculated from these combined trajectories correspond well to the experimental spectra for the buried site T4L F153R1 and the helix surface site T4L K65R1. This work provides a framework to further explore the modeling of the dynamics of the spin-label side chain in the wide variety of labeling environments encountered in site-directed spin labeling studies.     

Structure refinement of membrane proteins via molecular dynamics simulations

A refinement protocol based on physics‐based techniques established for water soluble proteins is tested for membrane protein structures. Initial structures were generated by homology modeling and sampled via molecular dynamics simulations in explicit lipid bilayer and aqueous solvent systems. Snapshots from the simulations were selected based on scoring with either knowledge‐based or implicit membrane‐based scoring functions and averaged to obtain refined models. The protocol resulted in consistent and significant refinement of the membrane protein structures similar to the performance of refinement methods for soluble proteins. Refinement success was similar between sampling in the presence of lipid bilayers and aqueous solvent but the presence of lipid bilayers may benefit the improvement of lipid‐facing residues. Scoring with knowledge‐based functions (DFIRE and RWplus) was found to be as good as scoring using implicit membrane‐based scoring functions suggesting that differences in internal packing is more important than orientations relative to the membrane during the refinement of membrane protein homology models.     

Lipid a from lipopolysaccharide recognition: Structure,dynamics and cooperativity by molecular dynamics simulations

Molecular dynamics simulations of Lipid A and its natural precursor Lipid IVA from E.coli have been carried out free in solution, bound to the myeliod differentiation protein 2 (MD2) and in the complex of MD2 with the toll like receptor 4 (TLR4). In addition, simulations of the ligand free MD2 and MD2‐TLR4 complex were performed. A structural and energetic characterization of the bound and unbound states of Lipid A/IVA was generated. As the crystal structures depict, the main driving force for MD2‐Lipid A/IVA are the hydrophobic interactions between the aliphatic tails and the MD2 cavity. The charged phosphate groups do strongly interact with positively charged residues, located at the surface of MD2. However, they are not essential for keeping the lipids in the cavity, indicating a more prominent role in binding recognition and ionic interactions with TLR4 at the MD2/TLR4 interface. Interestingly, in the absence of any ligand MD2 rapidly closes, blocking the binding cavity. The presence of TLR4, though changing the dynamics, was not able to impede the aforementioned closing event. We hypothesize that fluctuations of the H1 region are essential for this phenomenon, and it is plausible that an equilibrium between the open and closed states exists, although the lengths of our simulations are not sufficient to encompass the reversible process. The MD2/Lipid A‐TLR4 complex simulations show that the presence of the ligand energetically stabilizes the complex relative to the ligand‐free structures, indicating cooperativity in the binding process. © Proteins 2013. © 2012 Wiley Periodicals, Inc.     

Methods for molecular dynamics simulations of protein folding/unfolding in solution

All atom molecular dynamics simulations have become a standard method for mapping equilibrium protein dynamics and non-equilibrium events like folding and unfolding. Here, we present detailed methods for performing such simulations. Generic protocols for minimization, solvation, simulation, and analysis derived from previous studies are also presented. As a measure of validation, our water model is compared with experiment. An example of current applications of these methods, simulations of the ultrafast folding protein Engrailed Homeodomain are presented including the experimental evidence used to verify their results. Ultrafast folders are an invaluable tool for studying protein behavior as folding and unfolding events measured by experiment occur on timescales accessible with the high-resolution molecular dynamics methods we describe. Finally, to demonstrate the prospect of these methods for folding proteins, a temperature quench simulation of a thermal unfolding intermediate of the Engrailed Homeodomain is described.     

Structure, dynamics, and elasticity of free 16s rRNA helix 44 studied by molecular dynamics simulations

Molecular dynamics (MD) simulations were employed to investigate the structure, dynamics, and local base-pair step deformability of the free 16S ribosomal helix 44 from Thermus thermophilus and of a canonical A-RNA double helix. While helix 44 is bent in the crystal structure of the small ribosomal subunit, the simulated helix 44 is intrinsically straight. It shows, however, substantial instantaneous bends that are isotropic. The spontaneous motions seen in simulations achieve large degrees of bending seen in the X-ray structure and would be entirely sufficient to allow the dynamics of the upper part of helix 44 evidenced by cryo-electron microscopic studies. Analysis of local base-pair step deformability reveals a patch of flexible steps in the upper part of helix 44 and in the area proximal to the bulge bases, suggesting that the upper part of helix 44 has enhanced flexibility. The simulations identify two conformational substates of the second bulge area (bottom part of the helix) with distinct base pairing. In agreement with nuclear magnetic resonance (NMR) and X-ray studies, a flipped out conformational substate of conserved 1492A is seen in the first bulge area. Molecular dynamics (MD) simulations reveal a number of reversible alpha-gamma backbone flips that correspond to transitions between two known A-RNA backbone families. The flipped substates do not cumulate along the trajectory and lead to a modest transient reduction of helical twist with no significant influence on the overall geometry of the duplexes. Despite their considerable flexibility, the simulated structures are very stable with no indication of substantial force field inaccuracies.     

Structure of rhodamine 6G-DNA complexes from molecular dynamics simulations

Chromophore-DNA complexes are useful for understanding charge transport along pi-stacks once their structural properties have been clarified. We studied two rhodamine 6G semicapping complexes with 15-mer B-DNA duplexes to determine the preferred orientation of the dye with respect to the neighboring base pair. For each of these systems, two distinct chromophore alignments were identified and quantified in terms of base-step parameters. The obtained geometries agree well with those derived from an NMR structure refinement of similar complexes. Cross-correlation analysis of the base-step parameters shows that slide and twist are highly interdependent during the structural transition from one conformation to the other.     

Structure and dynamics of the full-length M1 muscarinic acetylcholine receptor studied by molecular dynamics simulations

The three-dimensional structure of full-length structure of the M₁ muscarinic receptor was obtained through the fragmental homology modeling procedure. A 10-ns molecular dynamics (MD) simulation of the protein imbedded in a lipid slab and surrounded by water molecules was further used to relax the model. It was found that the homology model corresponded to the conformation in the ground state, since no significant motions of the backbone of transmembrane domains were observed. Furthermore, the reliability of the model was validated by analyzing key inter-helical contacts, sidechain-sidechain interactions, the formation of stable aromatic microdomains (clusters) and the docking of acetylcholine to its binding site. Moreover, a few conserved interactions observed in the X-ray structure of rhodopsin, such as inter-helical sidechain-sidechain hydrogen bonds were accurately reproduced in the MD simulation. The coupling of ACh to its binding site was found to be dominated by π-cation and salt bridge interactions, while its conformational space was restrained through van der Waals and hydrogen bond interactions. In general, such features were in very good agreement with the available experimental as well as with theoretical data. Considering the above, the structural information obtained in this study can be used a starting point to investigate the activation mechanism of the receptor and the ability to develop selective agonists and allosteric modulators which could be used for the treatment of Alzheimer’s disease.     

Density functional and molecular dynamics simulations of local anesthetics in 0.9% NaCl solution

We have made density functional calculations and molecular dynamics (MD) simulations to investigate the structure and pharmacological action of local anesthetics: tetracaine, procaine and lidocaine. The MD simulations were made in a NPT ensemble, in a 0.9% NaCl solution, on both unprotonated and protonated forms of the molecules. The radial distribution function was used to study solvent effects in different regions of the molecules. Although all three anesthetics have different degrees of hydrophobicity, the amino-terminals were the mostly affected by the protonation yielding hydrophilic regions. The charged amino-esters present hydrophilicity on the ester as well as amine terminals. Cl^?  from the solvent solution forms hydrogen bonds via protonated hydrogen attached to nitrogen, yielding neutral molecules, which could, in principle, penetrate the membranes and loose Cl^?  to act in the protonated form. Density functional theory calculations indicated a change in the electrostatic potential and showed that Cl^?  weakly binds to the amine hydrogen, what suggests it is a favorable interaction and supports the existence of the hydrochloric forms of these local anesthetics.     

Structure and dynamics of interfacial water in an Lalpha phase lipid bilayer from molecular dynamics simulations

Based on molecular dynamics simulations, an analysis of structure and dynamics is performed on interfacial water at a liquid crystalline dipalmitoylphosphatidycholine/water system. Water properties relevant for understanding NMR relaxation are emphasized. The first and second rank orientational order parameters of the water O-H bonds were calculated, where the second rank order parameter is in agreement with experimental determined quadrupolar splittings. Also, two different interfacial water regions (bound water regions) are revealed with respect to different signs of the second rank order parameter. The water reorientation correlation function reveals a mixture of fast and slow decaying parts. The fast (ps) part of the correlation function is due to local anisotropic water reorientation whereas the much slower part is due to more complicated processes including lateral diffusion along the interface and chemical exchange between free and bound water molecules. The 100-ns-long molecular dynamics simulation at constant pressure (1 atm) and at a temperature of 50 degrees C of 64 lipid molecules and 64 x 23 water molecules lack a slow water reorientation correlation component in the ns time scale. The (2)H(2)O powder spectrum of the dipalmitoylphosphatidycholine/water system is narrow and consequently, the NMR relaxation time T(2) is too short compared to experimental results.     

Inhibitory mechanism of caspase-6 phosphorylation revealed by crystal structures, molecular dynamics simulations, and biochemical assays

The apoptotic effector caspase-6 (CASP6) has been clearly identified as a drug target due to its strong association with neurodegeneration and axonal pruning events as well as its crucial roles in Huntington disease and Alzheimer disease. CASP6 activity is suppressed by ARK5-mediated phosphorylation at Ser(257) with an unclear mechanism. In this work, we solved crystal structures of ΔproCASP6S257E and p20/p10S257E, which mimicked the phosphorylated CASP6 zymogen and activated CASP6, respectively. The structural investigation combined with extensive biochemical assay and molecular dynamics simulation studies revealed that phosphorylation on Ser(257) inhibited self-activation of CASP6 zymogen by "locking" the enzyme in the TEVD(193)-bound "inhibited state." The structural and biochemical results also showed that phosphorylation on Ser(257) inhibited the CASP6 activity by steric hindrance. These results disclosed the inhibition mechanism of CASP6 phosphorylation and laid the foundation for a new strategy of rational CASP6 drug design.     

Structure and interactions of RecA: plasticity revealed by molecular dynamics simulations

Eleven independent simulations, each involving three consecutive molecules in the RecA filament, carried out on the protein from Mycobacterium tuberculosis, Mycobacterium smegmatis and Escherichia coli and their Adenosine triphosphate (ATP) complexes, provide valuable information which is complementary to that obtained from crystal structures, in addition to confirming the robust common structural framework within which RecA molecules from different eubacteria function. Functionally important loops, which are largely disordered in crystal structures, appear to adopt in each simulation subsets of conformations from larger ensembles. The simulations indicate the possibility of additional interactions involving the P-loop which remains largely invariant. The phosphate tail of the ATP is firmly anchored on the loop while the nucleoside moiety exhibits substantial structural variability. The most important consequence of ATP binding is the movement of the ‘switch’ residue. The relevant simulations indicate the feasibility of a second nucleotide binding site, but the pathway between adjacent molecules in the filament involving the two nucleotide binding sites appears to be possible only in the mycobacterial proteins.     

Structure, interaction, dynamics and solvent effects on the DNA-EcoRI complex in aqueous solution from molecular dynamics simulation.

A 0.7-ns molecular dynamics simulation of the DNA-EcoRI complex in a 7.0-A solvent shell indicated a stable behavior of the system. No significant evaporation or smearing of the solvent's outer boundary occurred. The structure and the intermolecular interactions were found to be well maintained during the simulation. The interaction pattern in the simulation was found to be very similar to that in the crystal structure. Most of the specific interactions between the DNA and the protein were found to be enhanced in the simulation compared to that in the crystal structure as a result of improved interaction geometry. The nonspecific interactions were found to be stronger than the specific ones. The specific interactions between the N7 atoms of Gua(4) or Ade(5) or Ade(6) and the protein were found to be present over almost the entire time of the simulation, whereas hydrogen bonds involving the amino groups of the Ade(5) and Ade(6) with the protein were found to be relatively weaker, with lower probability and shorter lifetime. The time evolution of the root mean square deviations of the DNA and the protein were highly correlated even at the later part of the simulation, showing the tight binding between them. Several long-lived water bridges were found between the DNA backbone atoms and the protein and also between the two protein monomers, which increased the overall stability of the complex. The two protein monomers were found to interact strongly with each other. The energy of the DNA kink deformation was estimated as approximately 31 kcal/mol.