Unrestrained Molecular Dynamics Research Articles

Preventing lipophilic aggregation in cosolvent molecular dynamics simulations with hydrophobic probes using Plumed Automatic Restraining Tool (PART)

Cosolvent molecular dynamics (MD) simulations are molecular dynamics simulations used to identify preferable locations of small organic fragments on a protein target. Most cosolvent molecular dynamics workflows make use of only water-soluble fragments, as hydrophobic fragments would cause lipophilic aggregation. To date the two approaches that allow usage of hydrophobic cosolvent molecules are to use a low (0.2 M) concentration of hydrophobic probes, with the disadvantage of a lower sampling speed, or to use force field modifications, with the disadvantage of a difficult and inflexible setup procedure. Here we present a third alternative, that does not suffer from low sampling speed nor from cumbersome preparation procedures. We have built an easy-to-use open source command line tool PART (Plumed Automatic Restraining Tool) to generate a PLUMED file handling all intermolecular restraints to prevent lipophilic aggregation. We have compared restrained and unrestrained cosolvent MD simulations, showing that restraints are necessary to prevent lipophilic aggregation at hydrophobic probe concentrations of 0.5 M. Furthermore, we benchmarked PART generated restraints on a test set of four proteins (Factor-Xa, HIV protease, P38 MAP kinase and RNase A), showing that cosolvent MD with PART generated restraints qualitatively reproduces binding features of cocrystallised ligands.

Temperature-dependent elasticity of DNA, RNA, and hybrid double helices

Nucleic acid double helices in their DNA, RNA, and DNA-RNA hybrid form play a fundamental role in biology and are main building blocks of artificial nanostructures, but how their properties depend on temperature remains poorly understood. Here, we report thermal dependence of dynamic bending persistence length, twist rigidity, stretch modulus, and twist-stretch coupling for DNA, RNA, and hybrid duplexes between 7°C and 47°C. The results are based on all-atom molecular dynamics simulations using different force field parameterizations. We first demonstrate that unrestrained molecular dynamics can reproduce experimentally known mechanical properties of the duplexes at room temperature. Beyond experimentally known features, we also infer the twist rigidity and twist-stretch coupling of the hybrid duplex. As for the temperature dependence, we found that increasing temperature softens all the duplexes with respect to bending, twisting, and stretching. The relative decrease of the stretch moduli is 0.003–0.004/°C, similar for all the duplex variants despite their very different stretching stiffness, whereas RNA twist stiffness decreases by 0.003/°C, and smaller values are found for the other elastic moduli. The twist-stretch couplings are nearly unaffected by temperature. The stretching, bending, and twisting stiffness all include an important entropic component. Relation of our results to the two-state model of DNA flexibility is discussed. Our work provides temperature-dependent elasticity of nucleic acid duplexes at the microsecond scale relevant for initial stages of protein binding.

Exploring the disruption of SARS-CoV-2 RBD binding to hACE2.

The COVID-19 pandemic was declared due to the spread of the novel coronavirus, SARS-CoV-2. Viral infection is caused by the interaction between the SARS-CoV-2 receptor binding domain (RBD) and the human ACE2 receptor (hACE2). Previous computational studies have identified repurposed small molecules that target the RBD, but very few have screened drugs in the RBD-hACE2 interface. When studies focus solely on the binding affinity between the drug and the RBD, they ignore the effect of hACE2, resulting in an incomplete analysis. We screened ACE inhibitors and previously identified SARS-CoV-2 inhibitors for binding to the RBD-hACE2 interface, and then conducted 500ns of unrestrained molecular dynamics (MD) simulations of fosinopril, fosinoprilat, lisinopril, emodin, diquafosol, and physcion bound to the interface to assess the binding characteristics of these ligands. Based on MM-GBSA analysis, all six ligands bind favorably in the interface and inhibit the RBD-hACE2 interaction. However, when we repeat our simulation by first binding the drug to the RBD before interacting with hACE2, we find that fosinopril, fosinoprilat, and lisinopril result in a strongly interacting trimeric complex (RBD-drug-hACE2). Hydrogen bonding and pairwise decomposition analyses further suggest that fosinopril is the best RBD inhibitor. However, when lisinopril is bound, it stabilizes the trimeric complex and, therefore, is not an ideal potential drug candidate. Overall, these results reveal important atomistic interactions critical to the binding of the RBD to hACE2 and highlight the significance of including all protein partners in the evaluation of a potential drug candidate.

Molecular Motion in the Interconverting σ-H2, Di-, and Tri-hydride Regimes: .

The transition-metal complex is known to exist in three possible isomeric forms, including a nonclassical, σ-bound dihydrogen complex and two classical dihydride isomers. As such, it has served as a model complex for the energies of conversion between these limiting structural regimes. In the present study, ab initio molecular dynamics computer simulations, combined with enhanced sampling techniques, were utilized to directly assess the degree of motion and isomerization of the dihydrogen/dihydride moieties in this complex. Ligand rotations (for both the H2 unit and the phosphine units) were found to be dominant in the low-temperature (298 K) regime, and the classical thermodynamic distribution showed no probability of thermally accessing dihydride forms, although unrestrained molecular dynamics trajectories showed fleeting configurations outside of the σ-H2 configuration. Simulations at higher temperatures surprisingly revealed new tri-hydride isomers that are energetically competitive with the σ-H2 and cis-/trans-dihydride isomers. Low-energy pathways to hydrogen/hydride transfer and phosphine dissociation were readily accessible, which considerably expands the known isomeric flexibility of this complex.

Interactions between Rhodamine Dyes and Model Membrane Systems-Insights from Molecular Dynamics Simulations.

Background: rhodamines are dyes widely used as fluorescent tags in cell imaging, probing of mitochondrial membrane potential, and as P-glycoprotein model substrates. In all these applications, detailed understanding of the interaction between rhodamines and biomembranes is fundamental. Methods: we combined atomistic molecular dynamics (MD) simulations and fluorescence spectroscopy to characterize the interaction between rhodamines 123 and B (Rh123 and RhB, respectively) and POPC bilayers. Results: while the xanthene moiety orients roughly parallel to the membrane plane in unrestrained MD simulations, variations on the relative position of the benzoic ring (below the xanthene for Rh123, above it for RhB) were observed, and related to the structure of the two dyes and their interactions with water and lipids. Subtle distinctions were found among different ionization forms of the probes. Experimentally, RhB displayed a lipid/water partition coefficient more than two orders of magnitude higher than Rh123, in agreement with free energy profiles obtained from umbrella sampling MD. Conclusions: this work provided detailed insights on the similarities and differences in the behavior of bilayer-inserted Rh123 and RhB, related to the structure of the probes. The much higher affinity of RhB for the membranes increases the local concentration and explains its higher apparent affinity for P-glycoprotein reconstituted in model membranes.

Cryo-EM analysis provides new mechanistic insight into ATP binding to Ca2+ -ATPase SERCA2b.

Sarco/endoplasmic reticulum Ca2+ -ATPase (SERCA) 2b is a ubiquitous SERCA family member that conducts Ca2+ uptake from the cytosol to the ER. Herein, we present a 3.3Å resolution cryo-electron microscopy (cryo-EM) structure of human SERCA2b in the E1·2Ca2+ state, revealing a new conformation for Ca2+ -bound SERCA2b with a much closer arrangement of cytosolic domains than in the previously reported crystal structure of Ca2+ -bound SERCA1a. Multiple conformations generated by 3D classification of cryo-EM maps reflect the intrinsically dynamic nature of the cytosolic domains in this state. Notably, ATP binding residues of SERCA2b in the E1·2Ca2+ state are located at similar positions to those in the E1·2Ca2+ -ATP state; hence, the cryo-EM structure likely represents a preformed state immediately prior to ATP binding. Consistently, a SERCA2b mutant with an interdomain disulfide bridge that locks the closed cytosolic domain arrangement displayed significant autophosphorylation activity in the presence of Ca2+ . We propose a novel mechanism of ATP binding to SERCA2b.

Insight into the sequence-specific elements leading to increased DNA bending and ligase-mediated circularization propensity by antitumor trabectedin.

DNA curvature is the result of a combination of both intrinsic features of the double helix and external distortions introduced by the environment and the binding of proteins or drugs. The propensity of certain double-stranded DNA (dsDNA) sequences to bend is essential in crucial biological processes, such as replication and transcription, in which proteins are known to either recognize noncanonical DNA conformations or promote their formation upon DNA binding. Trabectedin (Yondelis®) is a clinically used antitumor drug which, following covalent bond formation with the 2-amino group of guanine, induces DNA curvature and enhances the circularization ratio, upon DNA ligation, of several dsDNA constructs but not others. By means of unrestrained molecular dynamics simulations using explicitly solvated all-atom models, we rationalize these experimental findings in structural terms and shed light on the crucial, albeit possibly underappreciated, role played by T4 DNA ligase in stabilizing a bent DNA conformation prior to cyclization. Taken together, our results expand our current understanding on how DNA shape modification by trabectedin may affect both the sequence-specific recognition by transcription factors to promoter sites and RNA polymerase II binding.

Structural Landscape of the Transition from an ssDNA Dumbbell Plus Its Complementary Hairpin to a dsDNA Microcircle Via a Kissing Loop Intermediate.

The experimental construction of a double-stranded DNA microcircle of only 42 base pairs entailed a great deal of ingenuity and hard work. However, figuring out the three-dimensional structures of intermediates and the final product can be particularly baffling. Using a combination of model building and unrestrained molecular dynamics simulations in explicit solvent we have characterized the different DNA structures involved along the process. Our 3D models of the single-stranded DNA molecules provide atomic insight into the recognition event that must take place for the DNA bases in the cohesive tail of the hairpin to pair with their complementary bases in the single-stranded loops of the dumbbell. We propose that a kissing loop involving six base pairs makes up the core of the nascent dsDNA microcircle. We also suggest a feasible pathway for the hybridization of the remaining complementary bases and characterize the final covalently closed dsDNA microcircle as possessing two well-defined U-turns. Additional models of the pre-ligation complex of T4 DNA ligase with the DNA dumbbell and the post-ligation pre-release complex involving the same enzyme and the covalently closed DNA microcircle are shown to be compatible with enzyme recognition and gap ligation.

Spontaneous rearrangement of acetylated xylan on hydrophilic cellulose surfaces

The interaction of xylan, an abundant plant polysaccharide, with cellulose microfibrils is essential for secondary cell wall strength. A deeper understanding of these interactions is crucial both to improve our understanding of plant cell wall architecture and to design alternate strategies to overcome cellulose recalcitrance for the production of biofuels and sustainable biomaterials. Naturally occurring acetate or glucuronic acid substitutions on xylan have been shown to influence xylan-cellulose interactions. Here, we use unrestrained molecular dynamics simulations to determine the interactions with the (110) hydrophilic face of cellulose fibers of four different xylans. In the absence of cellulose, all xylans, independent of the substitution pattern, adopt a highly flexible threefold helical screw conformation. However, when xylan is spatially close to a cellulose surface 1,2 linked acetyl xylans (2AcX) adopt rigid twofold helical screw conformations. The 2AcX conformations are primarily stabilized by interactions between the acetylated oxygen and the glycosidic linkage with C-O6 of cellulose. In contrast, the glycosidic oxygens and acetyl decorations for 1,3 linked acetyl groups (3AcX) are oriented away from the cellulose surface and the 3AcX xylans maintain threefold helical screw conformations on the cellulose surface. Our results show that evenly spaced chemical functionalization (with acetyl groups) and the position of substitution (1,2) on xylan backbone play key roles in tuning the xylan-cellulose interactions to stabilize the twofold helical screw conformations of xylan on the cellulose surface. A comparison with previous experimental findings further suggests that 1,2 substitutions induce twofold helical screw conformations of xylan on the cellulose surface irrespective of the chemical nature of the substituent, while 1,3 substitutions primarily bind lignin in threefold helical screw conformations rather than cellulose in plant cell walls.

Solution Structure of the Gloeobacter Violaceus Ligand-Gated Ion Channel Probed by Small-Angle Neutron Scattering

Pentameric ligand-gated ion channels transform chemical signals into electrical ones, a process during which they undergo subtle conformational cycling resulting in rapid, reversible gating of an intrinsic transmembrane pore. The pH-gated bacterial channel GLIC has proved a valuable model system for this receptor family, in part due to its accessibility to X-ray crystallography - which has provided structures of closed, open, and intermediate conformations. However, such models raise critical questions as to whether the X-ray structures represent actual functional states, and if the cycle includes additional conformations such as highly contracted or expanded states. Small-angle neutron scattering (SANS) offers an approach to characterizing solution-phase structure, increasingly applicable to biomolecular problems on this scale. Here we characterized GLIC by SANS, in combination with inline size-exclusion chromatography (SEC), match-out deuterated detergents, and computational modeling of multiple putative structural models. SEC-SANS under resting conditions (pH 7) enabled the characterization of non-aggegated species which fit predictions from resting-state models, particularly following unrestrained molecular dynamics simulation of X-ray structures. Curves collected under activating conditions (pH 3) best fit a mixture of resting- and open-state models, whereas they diverged from expanded-pore models. in addition to supporting less-expanded models of the open state, this work demonstrates that SEC-SANS, in combination with molecular modeling and simulations, can distinguish solution-phase functional states of an ion channel. As neutron source brilliance increases and inline SEC-SANS set-ups become increasingly accessible, such methods will offer valuable tools for the elucidation of receptor conformational cycling and pharmaceutical development.

Molecular dynamic simulations of diacridine binding to DNA: Indications that C6 diacridine can bisintercalate spanning two base pairs.

Dimers of 9-aminoacridine linked via the 9-amino group with polymethylene chains, termed diacridines, are known to bisintercalate into DNA when the linker comprises 6 or more methylene units. There are no literature reports of crystal or NMR solution structures for bisintercalated diacridine-DNA complexes, and the issue of the structure of the C6 ([CH2 ]n linker where n = 6) diacridine complex remains unresolved. Previously, based on simple geometric considerations, it was proposed that C6 diacridine could only span a single base pair, which requires that its bifunctional reaction violates the widely-observed "neighbor exclusion principle" where bound intercalators are separated by at least 2 base pairs. Here we have explored the structure of diacridine-DNA complexes using unrestrained molecular dynamics in explicit solvent using the parmbsc0 forcefield in AMBER14. We studied the C4 to C8 dimers, intercalated via both the minor and major DNA grooves, to a variety of nucleotide sequences. We find that C6, C7, and C8 diacridine are able to form 2 base pair bisintercalated complexes from either groove, whereas the C4 and C5 homologues cannot. We conclude that C6 diacridine does have the capacity to bisintercalate without violating neighbor exclusion, and that the previous proposed binding model needs revision.

On the role of solvent in hydrophobic cavity-ligand recognition kinetics.

A solvent often manifests itself as the key determinant of the kinetic aspect of the molecular recognition process. While the solvent is often depicted as a source of barrier in the ligand recognition process by the polar cavity, the nature of solvent's role in the recognition process involving hydrophobic cavity and hydrophobic ligand remains to be addressed. In this work, we quantitatively assess the role of solvent in dictating the kinetic process of recognition in a popular system involving the hydrophobic cavity and ligand. In this prototypical system, the hydrophobic cavity undergoes dewetting transition as the ligand approaches the cavity, which influences the cavity-ligand recognition kinetics. Here, we build a Markov state model (MSM) using adaptively sampled unrestrained molecular dynamics simulation trajectories to map the kinetic recognition process. The MSM-reconstructed free energy surface recovers a broad water distribution at an intermediate cavity-ligand separation, consistent with a previous report of dewetting transition in this system. Time-structured independent component analysis of the simulated trajectories quantitatively shows that cavity-solvent density contributes considerably in an optimized reaction coordinate involving cavity-ligand separation and water occupancy. Our approach quantifies two solvent-mediated macrostates at an intermediate separation of the cavity-ligand recognition pathways, apart from the fully ligand-bound and fully ligand-unbound macrostates. Interestingly, we find that these water-mediated intermediates, while transient in populations, can undergo slow mutual interconversion and create possibilities of multiple pathways of cavity recognition by the ligand. Overall, the work provides a quantitative assessment of the role that the solvent plays in facilitating the recognition process involving the hydrophobic cavity.

Permeation of nanoparticles across the intestinal lipid membrane: dependence on shape and surface chemistry studied through molecular simulations.

Nanoparticles are being explored for topical and oral drug delivery applications as they can cross various biological barriers, for example, the intestinal epithelium. The ability of nanoparticles to cross barriers depends on their morphological and surface properties such as size, surface chemistry and shape, among others. The effect of nanoparticle size on their membrane permeability has been well studied both experimentally and theoretically. However, less attention has been given to understand the role of nanoparticle shape in their translocation across biological barrier membranes. Here, we report on the influence of the nanoparticle's shape, surface chemistry and concentration on their permeation across a human intestinal apical cell membrane model. A representative multicomponent lipid bilayer model of the human intestinal apical membrane was built. The free energy of permeation of nanoparticles across the model lipid bilayer was calculated using multiple umbrella sampling simulations. The interaction of these nanoparticles with the model lipid bilayer was captured using extensive microsecond unrestrained molecular dynamics simulations. We observed that: (a) irrespective of the surface chemistry, the efficacy of nanoparticle penetration across the lipid layer was in the order of rod > disc > sphere; (b) irrespective of the shape, apolar and nonpolar nanoparticles were found to locate in the interior of the lipid bilayer, whereas charged and polar nanoparticles were either adsorbed on the lipid headgroups or remained in the water layer; (c) apolar and nonpolar disc shaped nanoparticles had higher efficacy in permeation across the lipid bilayer as compared to disc and sphere shaped nanoparticles; and (d) at a higher concentration of nanoparticles, sphere and disc shaped nanoparticles exhibited more agglomeration as compared to rod shaped nanoparticles. Based on these outcomes, a few nanoparticles were designed which penetrated readily into the lipid layer and these nanoparticles were also able to co-deliver a therapeutic protein inside the lipid layer. The apical model lipid membrane and protocols used in this study can thus be utilized for the in silico design of nanoparticles for the oral delivery of therapeutics.

Experimental and computational analysis of biased agonism on full-length and a C-terminally truncated adenosine A2A receptor

Biased agonism, the ability of agonists to differentially activate downstream signaling pathways by stabilizing specific receptor conformations, is a key issue for G protein-coupled receptor (GPCR) signaling. The C-terminal domain might influence this functional selectivity of GPCRs as it engages G proteins, GPCR kinases, β-arrestins, and several other proteins. Thus, the aim of this paper is to compare the agonist-dependent selectivity for intracellular pathways in a heterologous system expressing the full-length (A2AR) and a C-tail truncated (A2AΔ40R lacking the last 40 amino acids) adenosine A2A receptor, a GPCR that is already targeted in Parkinson’s disease using a first-in-class drug. Experimental data such as ligand binding, cAMP production, β-arrestin recruitment, ERK1/2 phosphorylation and dynamic mass redistribution assays, which correspond to different aspects of signal transduction, were measured upon the action of structurally diverse compounds (the agonists adenosine, NECA, CGS-21680, PSB-0777 and LUF-5834 and the SCH-58261 antagonist) in cells expressing A2AR and A2AΔ40R. The results show that taking cAMP levels and the endogenous adenosine agonist as references, the main difference in bias was obtained with PSB-0777 and LUF-5834. The C-terminus is dispensable for both G-protein and β-arrestin recruitment and also for MAPK activation. Unrestrained molecular dynamics simulations, at the μs timescale, were used to understand the structural arrangements of the binding cavity, triggered by these chemically different agonists, facilitating G protein binding with different efficacy.

Molecular mechanism of the skin permeation enhancing effect of ethanol: a molecular dynamics study

Ethanol is widely used in various pharmaceutical and cosmetic formulations in order to enhance skin penetration of active ingredients. While it is well known that ethanol partitions into the skin and enhances the permeation of both polar and nonpolar molecules, the exact mechanisms by which it enhances skin permeability are not fully understood. Several mechanisms have been proposed including lipid extraction from the stratum corneum (SC), fluidisation of SC lipid bilayer, alteration of SC protein conformation and enhancement of the drug solubility in the SC lipids. In this study, we performed molecular dynamics (MD) simulations of SC lipid bilayers comprised of an equimolar mixture of ceramides, cholesterol and free fatty acid in the presence of aqueous mixtures of ethanol. Various unrestrained MD simulations were performed in the presence of aqueous ethanol solution at molar ratios (x) ranging from x = 0 to x = 1. It was found that ethanol enhances bilayer permeability by dual actions (a) extraction of the skin lipids and (b) enhancing the mobility of lipid chains. Ethanol's permeation enhancing effect arises from its superior ability to form hydrogen bonds with headgroup atoms of skin lipids. Further, the free energy of extraction of ceramides (CER) and fatty acids (FFA) from the lipid bilayer was studied using umbrella sampling simulations. The free energy of extraction of CER was found to be much higher compared to FFA for all ethanol concentrations which shows that CER are difficult to extract as compared to FFA. Finally, the permeation of benzoic acid drug molecules through the skin lipid bilayer is shown in presence of ethanol molecules. It was found that ethanol selectively targets the FFA of the skin lipid bilayer and extracts it out of the lipid bilayer within few microseconds. Further, ethanol penetrates inside the lipid layer and creates the channels from which drug molecules can easily cross the lipid layer. Our observations (both in unrestrained and umbrella sampling simulations) are consistent with the experimental findings reported in the literature. The simulation methodology could be used for design and testing of permeation enhancers (acting on skin SC lipid lamella) for topical and transdermal drug delivery applications.

Crystallographic structure versus homology model: a case study of molecular dynamics simulation of human and zebrafish histone deacetylase 10

Histone deacetylase (HDAC) 10 has been implicated in the pathology of various cancers and neurodegenerative disorders, making the discovery of novel inhibitors of the isoform an important endeavor. However, the unavailability of crystallographic structure of human HDAC10 (hHDAC10) hinders structure-based drug design effort. Previously, we reported the homology modeled structure of human HDAC10 built using the crystallographic structure of Danio rerio (zebrafish) HDAC10 (zHDAC10) (Protein Data Bank (PDB) ID; 5TD7, released on 24 May 2017) as a template. Here, in continuation with our study, both hHDAC10 and zHDAC10, and their respective complexes with trichostatin A (TSA), quisinostat, and the native ligand (in 5TD7), 7-[(3-aminopropyl)amino]-1,1,1-trifluoroheptane-2,2-diol (PDB ID; FKS) were submitted to 100 ns-long unrestrained molecular dynamics (MD) simulations. Comparative analyses of the MD trajectories revealed that zHDAC10 and its complexes displayed higher stability than hHDAC10 and its corresponding complexes over time. Nonetheless, docking of active and inactive set molecules revealed that more reliable conformations of hHDAC10 could be obtained at an extended time period. This study may shed more light on the reliability of hHDAC10 modeled structure for use in selective inhibitor design.Communicated by Ramaswamy H. Sarma

Exploring Reactive Conformations of Coenzyme A during Binding and Unbinding to Pyruvate Formate-Lyase.

Pyruvate formate-lyase (PFL) is a glycyl radical enzyme that converts pyruvate and coenzyme A (CoA) into formate and acetyl-CoA in two half-reactions. Recently, we showed that the acetylation of the PFL active site in the first half-reaction induces subtle conformational changes, leading to the opening of a potential channel for CoA entry. Entry of CoA into the active site is crucial for the second half-reaction, involving the acetyl transfer to CoA, and the completion of the catalytic cycle. Using steered molecular dynamics (SMD) simulations, performed on acetylated and nonacetylated monomeric PFL model systems, we first of all investigate the possible entry/exit pathways of CoA with respect to the active site through the previously identified channel. We then perform umbrella sampling simulations on multiple snapshots from SMD trajectories as well as unrestrained molecular dynamics simulations starting from the final structures obtained from entry SMD, with a view to identifying possible bound states of CoA in the near vicinity of the active site. Detailed study of the unrestrained dissociation processes reveals the presence of stable and reactive bound states of CoA close to the active site, one of which is in an ideal position for triggering the second half-reaction. Examination of the spatial distributions associated with the reactive bound states allows us to discuss the free energy barriers. Umbrella sampling, performed on snapshots from unrestrained dynamics confirms the above findings. The significance of the results for the catalysis are discussed for both acetylated and nonacetylated systems.

Accuracy of MD solvent models in RNA structure refinement assessed via liquid-crystal NMR and spin relaxation data

Molecular dynamics (MD) simulations play an important role in characterizing Ribonucleic Acid (RNA) structure, augmenting information from experimental techniques such as Nuclear Magnetic Resonance (NMR). In this work, we examine the accuracy of structural representation resulting from application of a number of explicit and implicit solvent models and refinement protocols against experimental data ranging from high density of residual dipolar coupling (RDC) restraints to completely unrestrained simulations. For a prototype A-form RNA helix, our results indicate that AMBER RNA force field with either implicit or explicit solvent can produce a realistic dynamic representation of RNA helical structure, accurately cross-validating with respect to a diverse array of NMR observables. In refinement against NMR distance restraints, modern MD force fields are found to be equally adequate, with high fidelity cross-validation to the residual dipolar couplings (RDCs) and residual chemical shift anisotropies (RCSAs), while slightly over-estimating structural order as monitored via NMR relaxation data. With restraints trimmed to encode only for base pairing information, cross-validation quality significantly deteriorates, now exhibiting a pronounced dependence on the choice of the solvent model. This deterioration is found to be partially reversible by increasing planarity restraints on the nucleobase geometry. For completely unrestrained MD simulations, the choice of water model becomes very important, with the best-performing TIP4P-Ew accurately reproducing both the RDC and RCSA data, while closely matching the NMR-derived order parameters. The information provided here will serve as a foundation for MD-based refinement of solution state NMR structures of RNA.

A Molecular Dynamics Investigation of the Thermostability of Cold-Sensitive I707L KlenTaq1 DNA Polymerase and Its Wild-Type Counterpart.

DNA polymerase I from Thermus aquaticus ( Taq DNA polymerase) is useful for polymerase chain reactions because of its exceptional thermostability; however, its activity at low temperatures can cause amplification of unintended products. Mutation of isoleucine 707 to leucine (I707L) slows Taq DNA polymerase at low temperatures, which decreases unwanted amplification due to mispriming. In this work, unrestrained molecular dynamics (MD) simulations were performed on I707L and wild-type (WT) Taq DNA polymerase at 341 and 298 K to determine how the mutation affects the dynamic nature of the protein. The results suggest that I707L Taq DNA polymerase remains relatively immobile at room temperature and becomes more flexible at the higher temperature, while the WT Taq DNA polymerase demonstrates less substantial differences in dynamics at high and low temperatures. These results are in agreement with previous experimental results on the I707L mutant Taq DNA polymerase that show dynamic differences at high and low temperatures. The decreased mobility of the mutant at low temperature suggests that the mutant remains longer in the blocked conformation, and this may lead to reduced activity relative to the WT at 298 K. Principal component analysis revealed that the mutation results in decoupled movements of the Q helix and fingers domain. This decoupled nature of the mutant gives way to an increasingly flexible N-terminal end of the Q helix at 341 K, a characteristic not seen for WT Taq DNA polymerase.

Biomolecular Solvation Structure Revealed by Molecular Dynamics Simulations.

To compare ordered water positions from experiment with those from molecular dynamics (MD) simulations, a number of MD models of water structure in crystalline endoglucanase were calculated. The starting MD model was derived from a joint X-ray and neutron diffraction crystal structure, enabling the use of experimentally assigned protonation states. Simulations were performed in the crystalline state, using a periodic 2 × 2 × 2 supercell with explicit solvent. Water X-ray and neutron scattering density maps were computed from MD trajectories using standard macromolecular crystallography methods. In one set of simulations, harmonic restraints were applied to bias the protein structure toward the crystal structure. For these simulations, the recall of crystallographic waters using strong peaks in the MD water electron density was very good, and there also was substantial visual agreement between the boomerang-like wings of the neutron scattering density and the crystalline water hydrogen positions. An unrestrained simulation also was performed. For this simulation, the recall of crystallographic waters was much lower. For both restrained and unrestrained simulations, the strongest water density peaks were associated with crystallographic waters. The results demonstrate that it is now possible to recover crystallographic water structure using restrained MD simulations but that it is not yet reasonable to expect unrestrained MD simulations to do the same. Further development and generalization of MD water models for force-field development, macromolecular crystallography, and medicinal chemistry applications is now warranted. In particular, the combination of room-temperature crystallography, neutron diffraction, and crystalline MD simulations promises to substantially advance modeling of biomolecular solvation.