Sensitivity of molecular dynamics simulations to the choice of the X-ray structure used to model an enzymatic reaction

A subject of great practical importance that has not received much attention is the question of the sensitivity of molecular dynamics simulations to the initial X-ray structure used to set up the calculation. We have found two cases in which seemingly similar structures lead to quite different results, and in this article we present a detailed analysis of these cases. The first case is acyl-CoA dehydrogenase, and the chief difference of the two structures is attributed to a slight shift in a backbone carbonyl that causes a key residue (the proton-abstracting base) to be in a bad conformation for reaction. The second case is xylose isomerase, and the chief difference of the two structures appears to be the ligand sphere of a Mg2+ metal cofactor that plays an active role in catalysis.     

Conformational states of the glucocorticoid receptor DNA-binding domain from molecular dynamics simulations

Molecular dynamics simulations (MD) have been performed on variant crystal and NMR-derived structures of the glucocorticoid receptor DNA-binding domain (GR DBD). A loop region five residues long, the so-called D-box, exhibits significant flexibility, and transient perturbations of the tetrahedral geometry of two structurally important Cys4 zinc finger are seen, coupled to conformational changes in the D-box. In some cases, one of the Cys ligands to zinc exchanges with water, although no global distortion of the protein structure is observed. Thus, from MD simulation, dynamics of the D-box could partly be explained by solvent effects in conjunction with structural reformation of the zinc finger.     

Aggregation of small peptides studied by molecular dynamics simulations

Peptides and proteins tend to aggregate under appropriate conditions. The amyloid fibrils that are ubiquitously found among these structures are associated with major human diseases like Alzheimer's disease, type II diabetes, and various prion diseases. Lately, it has been observed that even very short peptides like tetra and pentapeptides can form ordered amyloid structures. Here, we present aggregation studies of three such small polypeptide systems, namely, the two amyloidogenic peptides DFNKF and FF, and a control (nonamyloidogenic) one, the AGAIL. The respective aggregation process is studied by all-atom Molecular Dynamics simulations, which allow to shed light on the fine details of the association and aggregation process. Our analysis suggests that naturally aggregating systems exhibit significantly diverse overall cluster shape properties and specific intermolecular interactions. Additional analysis was also performed on the previously studied NFGAIL system.     

Dynamics of the Hck-SH3 domain: comparison of experiment with multiple molecular dynamics simulations

Molecular dynamics calculations provide a method by which the dynamic properties of molecules can be explored over timescales and at a level of detail that cannot be obtained experimentally from NMR or X-ray analyses. Recent work (Philippopoulos M, Mandel AM, Palmer AG III, Lim C, 1997, Proteins 28:481-493) has indicated that the accuracy of these simulations is high, as measured by the correspondence of parameters extracted from these calculations to those determined through experimental means. Here, we investigate the dynamic behavior of the Src homology 3 (SH3) domain of hematopoietic cell kinase (Hck) via 5N backbone relaxation NMR studies and a set of four independent 4 ns solvated molecular dynamics calculations. We also find that molecular dynamics simulations accurately reproduce fast motion dynamics as estimated from generalized order parameter (S2) analysis for regions of the protein that have experimentally well-defined coordinates (i.e., stable secondary structural elements). However, for regions where the coordinates are not well defined, as indicated by high local root-mean-square deviations among NMR-determined structural family members or high B-factors/low electron density in X-ray crystallography determined structures, the parameters calculated from a short to moderate length (less than 5-10 ns) molecular dynamics trajectory are dependent on the particular coordinates chosen as a starting point for the simulation.     

Flexibility and charge asymmetry in the activation loop of Src tyrosine kinases

Regulated activity of Src kinases is critical for cell growth. Src kinases can be activated by trans-phosphorylation of a tyrosine located in the central activation loop of the catalytic domain. However, because the required exposure of this tyrosine is not observed in the down-regulated X-ray structures of Src kinases, transient partial opening of the activation loop appears to be necessary for such processes. Umbrella sampling molecular dynamics simulations are used to characterize the free energy landscape of opening of the hydrophilic part of the activation loop in the Src kinase Hck. The loop prefers a partially open conformation where Tyr416 has increased accessibility, but remains partly shielded. An asymmetric distribution of the charged residues in the sequence near Tyr416, which contributes to shielding, is found to be conserved in Src family members. A conformational equilibrium involving exchange of electrostatic interactions between the conserved residues Glu310 and Arg385 or Arg409 affects activation loop opening. A mechanism for access of unphosphorylated Tyr416 into an external catalytic site is suggested based on these observations.     

In silico studies of conformational dynamics of Mu opioid receptor performed using gaussian accelerated molecular dynamics

G-protein-coupled receptors play a crucial role in various signaling pathways and function as targets for treating a wide spectrum of diseases. Since the twentieth century, extensive research has been conducted on the Mu opioid receptor (MOR) as a drug target. We examined the MOR inactivation and activation processes using an enhanced sampling method (Gaussian accelerated molecular dynamics), the binding pocket site area method, the root mean square deviation method, and the free energy (potential of mean force) method. This study revealed two important intermediate MOR structures (intermediate and intermediate inactive), and the results suggest that the intermediate MOR structure is responsible for the selectivity of opioids.     

The active-inactive transition of human thymidylate synthase: targeted molecular dynamics simulations

Human thymidylate synthase (hTS) is an established anticancer target. It catalyses the production of 2'-deoxythymidine-5'-monophosphate, an essential building block for DNA synthesis. Because of the development of cellular drug resistance against current hTS inhibitors, alternative inhibition strategies are needed. hTS exists in two forms, active and inactive, defined by the conformation of the active-site (AS) loop, which carries the catalytic cysteine, C195. To investigate the mechanism of activation and inactivation, targeted molecular dynamics (TMD) simulations of the transitions between active and inactive states of hTS were performed. Analysis of changes in the dihedral angles in the AS loop during different TMD simulations revealed complex conformational transitions. Despite hTS being a homodimeric enzyme and the conformational transition significantly involving the dimer interface, the transition occurs in an asymmetric, sequential manner via an ensemble of pathways. In addition to C195, which reoriented during the simulations, other key residues in the rotation of the AS loop included W182 and R185. The interactions of the cognate bulky W182 residues at the dimer interface hindered the simultaneous twist of the AS loops in the hTS dimer. Interactions of R185, which is unique for hTS, with ligands at different allosteric sites affected the activation transition. In addition to providing insights into the activation/inactivation mechanism of hTS and how conformational transitions can occur in homodimeric proteins, our observations suggest that blocking of AS loop rotation by ligands binding in the large cavity between the loops could be one way to stabilize inactive hTS and inhibit the enzyme.     

Free energy and structural pathways of base flipping in a DNA GCGC containing sequence

Structural distortions of DNA are essential for its biological function due to the genetic information of DNA not being physically accessible in the duplex state. Base flipping is one of the simplest structural distortions of DNA and may represent an initial event in strand separation required to access the genetic code. Flipping is also utilized by DNA-modifying and repair enzymes to access specific bases. It is typically thought that base flipping (or base-pair opening) occurs via the major groove whereas minor groove flipping is only possible when mediated by DNA-binding proteins. Here, umbrella sampling with a novel center-of-mass pseudodihedral reaction coordinate was used to calculate the individual potentials of mean force (PMF) for flipping of the Watson-Crick (WC) paired C and G bases in the CCATGCGCTGAC DNA dodecamer. The novel reaction coordinate allowed explicit investigation of the complete flipping process via both the minor and major groove pathways. The minor and major groove barriers to flipping are similar for C base flipping while the major groove barrier is slightly lower for G base flipping. Minor groove flipping requires distortion of the WC partner while the flipping base pulls away from its partner during major groove flipping. The flipped states are represented by relatively flat free energy surfaces, with a small, local minimum observed for the flipped G base. Conserved patterns of phosphodiester backbone dihedral distortions during flipping indicate their essential role in the flipping process. During flipping, the target base tracks along the respective grooves, leading to hydrogen-bonding interactions with neighboring base-pairs. Such hydrogen-bonding interactions with the neighboring sequence suggest a novel mechanism of sequence dependence in DNA dynamics.     

Simulation of protein conformational freedom as a function of pH: constant-pH molecular dynamics using implicit titration

Solution pH is a determinant parameter on protein function and stability, and its inclusion in molecular dynamics simulations is attractive for studies at the molecular level. Current molecular dynamics simulations can consider pH only in a very limited way, through a somewhat arbitrary choice of a set of fixed charges on the titrable sites. Conversely, continuum electrostatic methods that explicitly treat pH effects assume a single protein conformation whose choice is not clearly defined. In this paper we describe a general method that combines both titration and conformational freedom. The method is based on a potential of mean force for implicit titration and combines both usual molecular dynamics and pH-dependent calculations based on continuum methods. A simple implementation of the method, using a mean field approximation, is presented and applied to the bovine pancreatic trypsin inhibitor. We believe that this constant-pH molecular dynamics method, by correctly sampling both charges and conformation, can become a valuable help in the understanding of the dependence of protein function and stability on pH. © 1997 Wiley-Liss Inc.     

Motional timescale predictions by molecular dynamics simulations: Case study using proline and hydroxyproline sidechain dynamics

We propose a new approach for force field optimizations which aims at reproducing dynamics characteristics using biomolecular MD simulations, in addition to improved prediction of motionally averaged structural properties available from experiment. As the source of experimental data for dynamics fittings, we use ¹³C NMR spin‐lattice relaxation times T₁ of backbone and sidechain carbons, which allow to determine correlation times of both overall molecular and intramolecular motions. For structural fittings, we use motionally averaged experimental values of NMR J couplings. The proline residue and its derivative 4‐hydroxyproline with relatively simple cyclic structure and sidechain dynamics were chosen for the assessment of the new approach in this work. Initially, grid search and simplexed MD simulations identified large number of parameter sets which fit equally well experimental J couplings. Using the Arrhenius‐type relationship between the force constant and the correlation time, the available MD data for a series of parameter sets were analyzed to predict the value of the force constant that best reproduces experimental timescale of the sidechain dynamics. Verification of the new force‐field (termed as AMBER99SB‐ILDNP) against NMR J couplings and correlation times showed consistent and significant improvements compared to the original force field in reproducing both structural and dynamics properties. The results suggest that matching experimental timescales of motions together with motionally averaged characteristics is the valid approach for force field parameter optimization. Such a comprehensive approach is not restricted to cyclic residues and can be extended to other amino acid residues, as well as to the backbone. Proteins 2014; 82:195–215. © 2013 Wiley Periodicals, Inc.     

The utility of molecular dynamics simulations for understanding site-directed mutagenesis of glycine residues in biotin carboxylase

Biotin carboxylase from Escherichia coli catalyzes the ATP-dependent carboxylation of biotin and is one component of the multienzyme complex acetyl-CoA carboxylase, which catalyzes the committed step in long-chain fatty acid synthesis. Comparison of the crystal structures of biotin carboxylase in the absence and presence of ATP showed a central B-domain closure when ATP was bound. Peptidic NH groups from two active site glycine residues (Gly165 and Gly166) that form hydrogen bonds to the phosphate oxygens of ATP were postulated to act as a "trigger" for movement of the B-domain. The function of these two glycine residues in the catalytic mechanism was studied by disruption of the hydrogen bonds using site-directed mutagenesis. Both single (G165V) and (G166V) and double mutants (G165V-G166V) were constructed. The mutations did not affect the maximal velocity of a partial reaction, the bicarbonate-dependent ATPase activity. This suggests that the peptidic NH groups of Gly165 and Gly166 are not triggers for domain movement. However, the K(m) values for ATP for each of the mutants was increased over 40-fold when compared with wild-type indicating the peptidic NH groups of Gly165 and Gly166 play a role in binding ATP. Consistent with ATP binding, the maximal velocity for the biotin-dependent ATPase activity (i.e. the complete reaction) was decreased over 100-fold suggesting the mutations have misaligned the reactants for optimal catalysis. Molecular dynamics studies confirm perturbation of the hydrogen bonds from the mutated residues to ATP, whereas the double mutant exhibits antagonistic effects such that hydrogen bonding from residues 165 and 166 to ATP is similar to that in the wild-type. Consistent with the site-directed mutagenesis results the molecular dynamics studies show that ATP is misaligned in the mutants.     

Folding simulations of Trp‐cage mini protein in explicit solvent using biasing potential replica‐exchange molecular dynamics simulations

Replica exchange molecular dynamics (RexMD) simulations are frequently used for studying structure formation and dynamics of peptides and proteins. A significant drawback of standard temperature RexMD is, however, the rapid increase of the replica number with increasing system size to cover a desired temperature range. A recently developed Hamiltonian RexMD method has been used to study folding of the Trp‐cage protein. It employs a biasing potential that lowers the backbone dihedral barriers and promotes peptide backbone transitions along the replica coordinate. In two independent applications of the biasing potential RexMD method including explicit solvent and starting from a completely unfolded structure the formation of near‐native conformations was observed after 30–40 ns simulation time. The conformation representing the most populated cluster at the final simulation stage had a backbone root mean square deviation of ～1.3 Å from the experimental structure. This was achieved with a very modest number of five replicas making it well suited for peptide and protein folding and refinement studies including explicit solvent. In contrast, during five independent continuous 70 ns molecular dynamics simulations formation of collapsed states but no near native structure formation was observed. The simulations predict a largely collapsed state with a significant helical propensity for the helical domain of the Trp‐cage protein already in the unfolded state. Hydrogen bonded bridging water molecules were identified that could play an active role by stabilizing the arrangement of the helical domain with respect to the rest of the chain already in intermediate states of the protein. Proteins 2009. © 2008 Wiley‐Liss, Inc.     

Effect of molecular confinement on internal enzyme dynamics: frequency domain fluorometry and molecular dynamics simulation studies

The tryptophanyl emission decay of the mesophilic beta-galactosidase from Aspergillus oryzae free in buffer and entrapped in agarose gel is investigated as a function of temperature and compared to that of the hyperthermophilic enzyme from Sulfolobus solfataricus. Both enzymes are tetrameric proteins with a large number of tryptophanyl residues, so the fluorescence emission can provide information on the conformational dynamics of the overall protein structure rather than that of the local environment. The tryptophanyl emission decays are best fitted by bimodal Lorentzian distributions. The long-lived component is ascribed to close, deeply buried tryptophanyl residues with reduced mobility; the short-lived one arises from tryptophanyl residues located in more flexible external regions of each subunit, some of which are involved in forming the catalytic site. The center of both lifetime distribution components at each temperature increases when going from the free in solution mesophilic enzyme to the gel-entrapped and hyperthermophilic enzyme, thus indicating that confinement of the mesophilic enzyme in the agarose gel limits the freedom of the polypeptide chain. A more complex dependence is observed for the distribution widths. Computer modeling techniques are used to recognize that the catalytic sites are similar for the mesophilic and hyperthermophilic beta-galactosidases. The effect due to gel entrapment is considered in dynamic simulations by imposing harmonic restraints to solvent-exposed atoms of the protein with the exclusion of those around the active site. The temperature dependence of the tryptophanyl fluorescence emission decay and the dynamic simulation confirm that more rigid structures, as in the case of the immobilized and/or hyperthermophilic enzyme, require higher temperatures to achieve the requisite conformational dynamics for an effective catalytic action and strongly suggest a link between conformational rigidity and enhanced thermal stability.     

Modeling of the structure in aqueous solution of the exopolysaccharide produced by Lactobacillus helveticus 766.

A method is described for constructing a conformational model in water of a heteropolysaccharide built up from repeating units, and is applied to the exopolysaccharide produced by Lactobacillus helveticus 766. The molecular modeling method is based on energy minima, obtained from molecular mechanics calculations of each of the constituting disaccharide fragments of the repeating unit in vacuo, as starting points. Subsequently, adaptive umbrella sampling of the potential of mean force is applied to extract rotamer populations of glycosidic dihedral angles of oligosaccharide fragments in solution. From these analyses, the most probable conformations are constructed for the hexasaccharide-repeating unit of the polysaccharide. After exploring the conformational space of each of the individual structures by molecular dynamics simulations, the different repeating unit conformations are used as building blocks for the generation of oligo- and polysaccharide models, by using a polysaccharide building program. The created models of the exopolysaccharide produced by L. helveticus 766 exhibit a flexible twisted secondary structure and tend to adopt a random coil conformation as tertiary structure.     

Possible ligand release pathway of dipeptidyl peptidase IV investigated by molecular dynamics simulations

Dipeptidyl peptidase IV (DPP4) is an important target for the treatment of Type II diabetes mellitus. The crystal structure of DPP4 demonstrates that there are two possible pathways to the active site, a side opening and a β propeller opening. However, it still lacks quantitative evidence to illustrate which pathway is more favorable for inhibitor to enter into or release from the active site. In this study, conventional and steered molecular dynamics simulations were performed to explore the details of inhibitor Q448 release from the active site of DPP4 via the two potential pathways. The comparisons of force and work together with potentials of mean force results suggested that the side opening might be more favorable for the inhibitor to pass through. Moreover, Glu205-Glu206 and Phe357 were recognized as two "key residues" in the active site for inhibitor binding. Accordingly, suggestions for further inhibitor design were provided.     

Conformation of the Ras-binding domain of Raf studied by molecular dynamics and free energy simulations

Recognition of Ras by its downstream target Raf is mediated by a Ras-recognition region in the Ras-binding domain (RBD) of Raf. Residues 78–89 in this region occupy two different conformations in the ensemble of NMR solution structures of the RBD: a fully α-helical one, and one where 87–90 form a type IV β-turn. Molecular dynamics simulations of the RBD in solution were performed to explore the stability of these and other possible conformations of both the wild-type RBD and the R89K mutant, which does not bind Ras. The simulations sample a fully helical conformation for residues 78–89 similar to the NMR helical structures, a conformation where 85–89 form a 3₁₀-helical turn, and a conformation where 87–90 form a type I |iB-turn, whose free energies are all within 0.3 kcal/mol of each other. NOE patterns and H^α chemical shifts from the simulations are in reasonable agreement with experiment. The NMR turn structure is calculated to be 3 kcal/mol higher than the three above conformations. In a simulation with the same implicit solvent model used in the NMR structure generation, the turn conformation relaxes into the fully helical conformation, illustrating possible structural artifacts introduced by the implicit solvent model. With the Raf R89K mutant, simulations sample a fully helical and a turn conformation, the turn being 0.9 kcal/mol more stable. Thus, the mutation affects the population of RBD conformations, and this is expected to affect Ras binding. For example, if the fully helical conformation of residues 78–89 is required for binding, its free energy increase in R89K will increase the binding free energy by about 0.6 kcal/mol. Proteins 31:186–200, 1998. © 1998 Wiley-Liss, Inc.     

Multiple conformational states and gate opening of outer membrane protein TolC revealed by molecular dynamics simulations

Outer membrane protein TolC serves as an exit duct for exporting substances out of cell. The occluded periplasmic entrance of TolC is required to open for substrate transport, although the opening mechanism remains elusive. In this study, systematic molecular dynamics (MD) simulations for wild type TolC and six mutants were performed to explore the conformational dynamics of TolC. The periplasmic gate was shown to sample multiple conformational states with various degrees of gating opening. The gate opening was facilitated by all mutations except Y362F, which adopts an even more closed state than wild type TolC. The interprotomer salt‐bridge R367–D153 is turned out to be crucial for periplasmic gate opening. The mutations that disrupt the interactions at the periplasmic tip may affect the stability of the trimeric assembly of TolC. Structural asymmetry of the periplasmic gate was observed to be opening size dependent. Asymmetric conformations are found in moderately opening states, while the most and the least opening states are often more symmetric. Finally, it is shown that lowering pH can remarkably stabilize the closed state of the periplasmic gate. Proteins 2014; 82:2169–2179. © 2014 Wiley Periodicals, Inc.