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.     

Improving the accuracy of protein stability predictions with multistate design using a variety of backbone ensembles

Multistate computational protein design (MSD) with backbone ensembles approximating conformational flexibility can predict higher quality sequences than single‐state design with a single fixed backbone. However, it is currently unclear what characteristics of backbone ensembles are required for the accurate prediction of protein sequence stability. In this study, we aimed to improve the accuracy of protein stability predictions made with MSD by using a variety of backbone ensembles to recapitulate the experimentally measured stability of 85 Streptococcal protein G domain β1 sequences. Ensembles tested here include an NMR ensemble as well as those generated by molecular dynamics (MD) simulations, by Backrub motions, and by PertMin, a new method that we developed involving the perturbation of atomic coordinates followed by energy minimization. MSD with the PertMin ensembles resulted in the most accurate predictions by providing the highest number of stable sequences in the top 25, and by correctly binning sequences as stable or unstable with the highest success rate (≈90%) and the lowest number of false positives. The performance of PertMin ensembles is due to the fact that their members closely resemble the input crystal structure and have low potential energy. Conversely, the NMR ensemble as well as those generated by MD simulations at 500 or 1000 K reduced prediction accuracy due to their low structural similarity to the crystal structure. The ensembles tested herein thus represent on‐ or off‐target models of the native protein fold and could be used in future studies to design for desired properties other than stability. Proteins 2014; 82:771–784. © 2013 Wiley Periodicals, Inc.     

Selective refinement and selection of near‐native models in protein structure prediction

In recent years in silico protein structure prediction reached a level where fully automated servers can generate large pools of near‐native structures. However, the identification and further refinement of the best structures from the pool of models remain problematic. To address these issues, we have developed (i) a target‐specific selective refinement (SR) protocol; and (ii) molecular dynamics (MD) simulation based ranking (SMDR) method. In SR the all‐atom refinement of structures is accomplished via the Rosetta Relax protocol, subject to specific constraints determined by the size and complexity of the target. The best‐refined models are selected with SMDR by testing their relative stability against gradual heating through all‐atom MD simulations. Through extensive testing we have found that Mufold‐MD, our fully automated protein structure prediction server updated with the SR and SMDR modules consistently outperformed its previous versions. Proteins 2015; 83:1823–1835. © 2015 Wiley Periodicals, Inc.     

Application of biasing‐potential replica‐exchange simulations for loop modeling and refinement of proteins in explicit solvent

Comparative or homology modeling of a target protein based on sequence similarity to a protein with known structure is widely used to provide structural models of proteins. Depending on the target‐template similarity these model structures may contain regions of limited structural accuracy. In principle, molecular dynamics (MD) simulations can be used to refine protein model structures and also to model loop regions that connect structurally conserved regions but it is limited by the currently accessible simulation time scales. A recently developed biasing potential replica exchange (BP‐REMD) method was used to refine loops and complete decoy protein structures at atomic resolution including explicit solvent. In standard REMD simulations several replicas of a system are run in parallel at different temperatures allowing exchanges at preset time intervals. In a BP‐REMD simulation replicas are controlled by various levels of a biasing potential to reduce the energy barriers associated with peptide backbone dihedral transitions. The method requires much fewer replicas for efficient sampling compared with T‐REMD. Application of the approach to several protein loops indicated improved conformational sampling of backbone dihedral angle of loop residues compared to conventional MD simulations. BP‐REMD refinement simulations on several test cases starting from decoy structures deviating significantly from the native structure resulted in final structures in much closer agreement with experiment compared to conventional MD simulations. Proteins 2010. © 2010 Wiley‐Liss, Inc.     

Restrained molecular dynamics of solvated duplex DNA using the particle mesh Ewald method

Restrained and unrestrained aqueous solution molecular dynamics simulations applying the particle mesh Ewald (PME) method to DNA duplex structures previously determined via in vacuo restrained molecular dynamics with NMR-derived restraints are reported. Without experimental restraints, the DNA decamer, d(CATTTGCATC)d(GATGCAAATG) and trisdecamer, d(AGCTTGCCTTGAG)d(CTCAAGGCAAGCT), structures are stable on the nanosecond time scale and adopt conformations in the B-DNA family. These free DNA simulations exhibit behavior characteristic of PME simulations previously performed on DNA sequences, including a low helical twist, frequent sugar pucker transitions, BI- BII(–) transitions and coupled crankshaft (–) motion. Refinement protocols similar to the original in vacuo restrained molecular dynamics (RMD) refinements but in aqueous solution using the Cornell et al. force field [Cornell et al. (1995) J. Am. Chem. Soc., 117, 5179–5197] and a particle mesh Ewald treatment produce structures which fit the restraints very well and are very similar to the original in vacuo NMR structure, except for a significant difference in the average helical twist. Figures of merit for the average structure found in the RMD PME decamer simulations in solution are equivalent to the original in vacuo NMR structure while the figures of merit for the free MD simulations are significantly higher. The free MD simulations with the PME method, however, lead to some sequence-dependent structural features in common with the NMR structures, unlike free MD calculations with earlier force fields and protocols. There is some suggestion that the improved handling of electrostatics by PME improves long-range structural aspects which are not well defined by the short-range nature of NMR restraints.     

In silico screening,molecular docking,and molecular dynamics studies of SNP-derived human P5CR mutants

Pyrroline-5-carboxylate reductase (P5CR) encoded by PYCR1 gene is a housekeeping enzyme that catalyzes the reduction of P5C to proline using NAD(P)H as the cofactor. In this study, we used in silico approaches to examine the role of nonsynonymous single-nucleotide polymorphisms in the PYCR1 gene and their putative functions in the pathogenesis of Cutis Laxa. Among the 348 identified SNPs, 15 were predicted to be potentially damaging by both SIFT and PolyPhen tools; of them two SNP‐derived mutations, R119G and G206W, have been previously reported to correlate with Cutis Laxa. These two mutations were therefore selected to be mapped to the wild‐type (WT) P5CR structure for further structural and functional analyses. The results of comparative computational analyses using I-Mutant and Autodock reveal reductions in both stability and cofactor binding affinity of these two mutants. Comparative molecular dynamics (MD) simulations were performed to evaluate the changes in dynamic properties of P5CR upon mutations. The results reveal that the two mutations enhance the rigidity of P5CR structure, especially that of cofactor binding site, which could result in decreased kinetics of cofactor entrance and egress. Comparison between the structural properties of the WT and mutants during MD simulations shows that the enhanced rigidity of mutants results most likely from the increased number of inter‐atomic interactions and the decreased number of dynamic hydrogen bonds. Our study provides novel insight into the deleterious effects of the R119G and G206W mutations on P5CR, and sheds light on the mechanisms by which these mutations mediate Cutis Laxa.     

Refinement of protein structure homology models via long, all-atom molecular dynamics simulations

Accurate computational prediction of protein structure represents a longstanding challenge in molecular biology and structure-based drug design. Although homology modeling techniques are widely used to produce low-resolution models, refining these models to high resolution has proven difficult. With long enough simulations and sufficiently accurate force fields, molecular dynamics (MD) simulations should in principle allow such refinement, but efforts to refine homology models using MD have for the most part yielded disappointing results. It has thus far been unclear whether MD-based refinement is limited primarily by accessible simulation timescales, force field accuracy, or both. Here, we examine MD as a technique for homology model refinement using all-atom simulations, each at least 100 μs long-more than 100 times longer than previous refinement simulations-and a physics-based force field that was recently shown to successfully fold a structurally diverse set of fast-folding proteins. In MD simulations of 24 proteins chosen from the refinement category of recent Critical Assessment of Structure Prediction (CASP) experiments, we find that in most cases, simulations initiated from homology models drift away from the native structure. Comparison with simulations initiated from the native structure suggests that force field accuracy is the primary factor limiting MD-based refinement. This problem can be mitigated to some extent by restricting sampling to the neighborhood of the initial model, leading to structural improvement that, while limited, is roughly comparable to the leading alternative methods.     

Structure–function relationships in ShKT domain peptides: ShKT-Ts1 from the sea anemone Telmatactis stephensoni

Diverse structural scaffolds have been described in peptides from sea anemones, with the ShKT domain being a common scaffold first identified in ShK toxin from Stichodactyla helianthus. ShK is a potent blocker of voltage-gated potassium channels (K_V1.x), and an analog, ShK-186 (dalazatide), has completed Phase 1 clinical trials in plaque psoriasis. The ShKT domain has been found in numerous other species, but only a tiny fraction of ShKT domains has been characterized functionally. Despite adopting the canonical ShK fold, some ShKT peptides from sea anemones inhibit K_V1.x, while others do not. Mutagenesis studies have shown that a Lys–Tyr (KY) dyad plays a key role in K_V1.x blockade, although a cationic residue followed by a hydrophobic residue may also suffice. Nevertheless, ShKT peptides displaying an ShK-like fold and containing a KY dyad do not necessarily block potassium channels, so additional criteria are needed to determine whether new ShKT peptides might show activity against potassium channels. In this study, we used a combination of NMR and molecular dynamics (MD) simulations to assess the potential activity of a new ShKT peptide. We determined the structure of ShKT-Ts1, from the sea anemone Telmatactis stephensoni, examined its tissue localization, and investigated its activity against a range of ion channels. As ShKT-Ts1 showed no activity against K_V1.x channels, we used MD simulations to investigate whether solvent exposure of the dyad residues may be informative in rationalizing and potentially predicting the ability of ShKT peptides to block K_V1.x channels. We show that either a buried dyad that does not become exposed during MD simulations, or a partially exposed dyad that becomes buried during MD simulations, correlates with weak or absent activity against K_V1.x channels. Therefore, structure determination coupled with MD simulations, may be used to predict whether new sequences belonging to the ShKT family may act as potassium channel blockers.     

A molecular dynamics study of a miRNA:mRNA interaction

In this paper we present a methodology to evaluate the binding free energy of a miRNA:mRNA complex through molecular dynamics (MD)–thermodynamic integration (TI) simulations. We applied our method to the Caenorhabditis elegans let-7 miRNA:lin-41 mRNA complex—a validated miRNA:mRNA interaction—in order to estimate the energetic stability of the structure. To make the miRNA:mRNA simulation possible and realistic, the methodology introduces specific solutions to overcome some of the general challenges of nucleic acid simulations and binding free energy computations that have been discussed widely in many previous research reports. The main features of the proposed methodology are: (1) positioning of the restraints imposed on the simulations in order to guarantee complex stability; (2) optimal sampling of the phase space to achieve satisfactory accuracy in the binding energy value; (3) determination of a suitable trade-off between computational costs and accuracy of binding free energy computation by the assessment of the scalability characteristics of the parallel simulations required for the TI. The experiments carried out demonstrate that MD simulations are a viable strategy for the study of miRNA binding characteristics, opening the way to the development of new computational target prediction methods based on three-dimensional structure information.     

Observation of an A-DNA to B-DNA transition in a nonhelical nucleic acid hairpin molecule using molecular dynamics.

One of the truly challenging problems for molecular dynamics (MD) simulations is demonstrating that the trajectories can sample not only in the vicinity of an experimentally determined structure, but also that the trajectories can find the correct experimental structure starting from some other structure. Frequently these transitions to the correct structure require that the simulations overcome energetic barriers to conformational change. Here we present unrestrained molecular dynamics simulations of the DNA analogs of the RNA 5'-GGACUUCGGUCC-3' hairpin tetraloop. In one simulation we have used deoxyuracil residues, and in the other we have used the native DNA deoxythymine residues. We demonstrate that, on a nanosecond time scale, MD is able to simulate the transitions of both of the A-DNA stems to B-DNA stems within the constraints imposed by the four-base loop that caps the helix. These results suggest that we are now in a position to use MD to address the nature of sequence-dependent structural effects in nonduplex DNA structures.     

A novel method for predicting and using distance constraints of high accuracy for refining protein structure prediction

The principal bottleneck in protein structure prediction is the refinement of models from lower accuracies to the resolution observed by experiment. We developed a novel constraints‐based refinement method that identifies a high number of accurate input constraints from initial models and rebuilds them using restrained torsion angle dynamics (rTAD). We previously created a Bayesian statistics‐based residue‐specific all‐atom probability discriminatory function (RAPDF) to discriminate native‐like models by measuring the probability of accuracy for atom type distances within a given model. Here, we exploit RAPDF to score (i.e., filter) constraints from initial predictions that may or may not be close to a native‐like state, obtain consensus of top scoring constraints amongst five initial models, and compile sets with no redundant residue pair constraints. We find that this method consistently produces a large and highly accurate set of distance constraints from which to build refinement models. We further optimize the balance between accuracy and coverage of constraints by producing multiple structure sets using different constraint distance cutoffs, and note that the cutoff governs spatially near versus distant effects in model generation. This complete procedure of deriving distance constraints for rTAD simulations improves the quality of initial predictions significantly in all cases evaluated by us. Our procedure represents a significant step in solving the protein structure prediction and refinement problem, by enabling the use of consensus constraints, RAPDF, and rTAD for protein structure modeling and refinement. Proteins 2009. © 2009 Wiley‐Liss, Inc.     

Interactions Between a Voltage Sensor and a Toxin via Multiscale Simulations

Gating-modifier toxins inhibit voltage-gated ion channels by binding the voltage sensors (VS) and altering the energetics of voltage-dependent gating. These toxins are thought to gain access to the VS via the membrane (i.e., by partitioning from water into the membrane before binding the VS). We used serial multiscale molecular-dynamics (MD) simulations, via a combination of coarse-grained (CG) and atomistic (AT) simulations, to study how the toxin VSTx1, which inhibits the archeabacterial voltage-gated potassium channel KvAP, interacts with an isolated membrane-embedded VS domain. In the CG simulations, VSTx1, which was initially located in water, partitioned into the headgroup/water interface of the lipid bilayer before binding the VS. The CG configurations were used to generate AT representations of the system, which were subjected to AT-MD to further evaluate the stability of the complex and refine the predicted VS/toxin interface. VSTx1 interacted with a binding site on the VS formed by the C-terminus of S1, the S1-S2 linker, and the N-terminus of S4. The predicted VS/toxin interactions are suggestive of toxin-mediated perturbations of the interaction between the VS and the pore domain of Kv channels, and of the membrane. Our simulations support a membrane-access mechanism of inhibition of Kv channels by VS toxins. Overall, the results show that serial multiscale MD simulations may be used to model a two-stage process of protein-bilayer and protein-protein interactions within a membrane.     

A molecular dynamics study of the C-terminal fragment of the L7/L12 ribosomal protein

The crystallographic dimer of the C-terminal fragment (CTF) of the L7/L12 ribosomal protein has been subjected to molecular dynamics (MD) simulations. A 90 picosecond (ps) trajectory for the protein dimer, 19 water molecules and two counter ions has been calculated at constant temperature. Effects of intermolecular interactions on the structure and dynamics have been studied. The exact crystallographic symmetry is lost and the atomic fluctuations differ from one monomer to the other. The average MD structure is more stable than the X-ray one, as judged by accessible surface area and energy calculations. Crystal (non-dimeric) interactions have been simulated in another 40 ps trajectory by using harmonic restraints to represent intermolecular hydrogen bonds. The conformational changes with respect ot the X-ray structure are then virtually suppressed.The unrestrained dimer trajectory has been scanned for cooperative motions involving secondary structure elements. The intrinsic collective motions of the monomer are transmitted via intermolecular contacts to the dimer structure.The existence of a stable dimeric form of CTF, resembling the crystallographic one, has been documented. At the cost of fairly small energy expenditure the dimer has considerable conformational flexibility. This flexibility may endow the dimer with some functional potential as an energy transducer.     

Control of cation permeation through the nicotinic receptor channel

We used molecular dynamics (MD) simulations to explore the transport of single cations through the channel of the muscle nicotinic acetylcholine receptor (nAChR). Four MD simulations of 16 ns were performed at physiological and hyperpolarized membrane potentials, with and without restraints of the structure, but all without bound agonist. With the structure unrestrained and a potential of −100 mV, one cation traversed the channel during a transient period of channel hydration; at −200 mV, the channel was continuously hydrated and two cations traversed the channel. With the structure restrained, however, cations did not traverse the channel at either membrane potential, even though the channel was continuously hydrated. The overall results show that cation selective transport through the nAChR channel is governed by electrostatic interactions to achieve charge selectivity, but ion translocation relies on channel hydration, facilitated by a trans-membrane field, coupled with dynamic fluctuations of the channel structure.     

Temperature‐accelerated molecular dynamics gives insights into globular conformations sampled in the free state of the AC catalytic domain

The catalytic domain of the adenyl cyclase (AC) toxin from Bordetella pertussis is activated by interaction with calmodulin (CaM), resulting in cAMP overproduction in the infected cell. In the X‐ray crystallographic structure of the complex between AC and the C terminal lobe of CaM, the toxin displays a markedly elongated shape. As for the structure of the isolated protein, experimental results support the hypothesis that more globular conformations are sampled, but information at atomic resolution is still lacking. Here, we use temperature‐accelerated molecular dynamics (TAMD) simulations to generate putative all‐atom models of globular conformations sampled by CaM‐free AC. As collective variables, we use centers of mass coordinates of groups of residues selected from the analysis of standard molecular dynamics (MD) simulations. Results show that TAMD allows extended conformational sampling and generates AC conformations that are more globular than in the complexed state. These structures are then refined via energy minimization and further unrestrained MD simulations to optimize inter‐domain packing interactions, thus resulting in the identification of a set of hydrogen bonds present in the globular conformations. Proteins 2014; 82:2483–2496. © 2014 Wiley Periodicals, Inc.     

The three‐dimensional structure of TrmB,a transcriptional regulator of dual function in the hyperthermophilic archaeon Pyrococcus furiosus in complex with sucrose

TrmB is a repressor that binds maltose, maltotriose, and sucrose, as well as other α‐glucosides. It recognizes two different operator sequences controlling the TM (Trehalose/Maltose) and the MD (Maltodextrin) operon encoding the respective ABC transporters and sugar‐degrading enzymes. Binding of maltose to TrmB abrogates repression of the TM operon but maintains the repression of the MD operon. On the other hand, binding of sucrose abrogates repression of the MD operon but maintains repression of the TM operon. The three‐dimensional structure of TrmB in complex with sucrose was solved and refined to a resolution of 3.0 Å. The structure shows the N‐terminal DNA binding domain containing a winged‐helix‐turn‐helix (wHTH) domain followed by an amphipathic helix with a coiled‐coil motif. The latter promotes dimerization and places the symmetry mates of the putative recognition helix in the wHTH motif about 30 Å apart suggesting a canonical binding to two successive major grooves of duplex palindromic DNA. This suggests that the structure resembles the conformation of TrmB recognizing the pseudopalindromic TM promoter but not the conformation recognizing the nonpalindromic MD promoter.     

Can molecular dynamics simulations improve the structural accuracy and virtual screening performance of GPCR models?

The determination of G protein-coupled receptor (GPCR) structures at atomic resolution has improved understanding of cellular signaling and will accelerate the development of new drug candidates. However, experimental structures still remain unavailable for a majority of the GPCR family. GPCR structures and their interactions with ligands can also be modelled computationally, but such predictions have limited accuracy. In this work, we explored if molecular dynamics (MD) simulations could be used to refine the accuracy of in silico models of receptor-ligand complexes that were submitted to a community-wide assessment of GPCR structure prediction (GPCR Dock). Two simulation protocols were used to refine 30 models of the D₃ dopamine receptor (D₃R) in complex with an antagonist. Close to 60 μs of simulation time was generated and the resulting MD refined models were compared to a D₃R crystal structure. In the MD simulations, the receptor models generally drifted further away from the crystal structure conformation. However, MD refinement was able to improve the accuracy of the ligand binding mode. The best refinement protocol improved agreement with the experimentally observed ligand binding mode for a majority of the models. Receptor structures with improved virtual screening performance, which was assessed by molecular docking of ligands and decoys, could also be identified among the MD refined models. Application of weak restraints to the transmembrane helixes in the MD simulations further improved predictions of the ligand binding mode and second extracellular loop. These results provide guidelines for application of MD refinement in prediction of GPCR-ligand complexes and directions for further method development.